ecearth3/sources/tm5mp/base/tmm_mf_tm5_nc.F90

!###############################################################################
!
! Input/output of meteofiles produced by TM5 : NetCDF 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_nc

  use GO, only : gol, goErr, goPr
  
  use GO , only : TDate
  !use MDF, only : MDF_NETCDF
  use MDF, only : MDF_NETCDF4

  implicit none
  

  ! --- in/out -----------------------------------
  
  private
  
  public  ::  TMeteoFile_tm5_nc
  public  ::  TMM_MF_TM5_NC_Init, TMM_MF_TM5_NC_Done
  public  ::  Init, Done
  public  ::  Get
  public  ::  ReadRecord
  public  ::  WriteRecord
  

  ! --- const ------------------------------------
  
  character(len=*), parameter  ::  mname = 'tmm_mf_tm5_nc'
  
  ! ~~~ output keys and defaults
  
  ! current format version
  character(len=*), parameter  ::  output_format = 'tm5-nc'
  
  ! extension and type:
  !integer, parameter           ::  output_type = MDF_NETCDF
  integer, parameter           ::  output_type = MDF_NETCDF4
  character(len=*), parameter  ::  output_ext  = '.nc'
  
  ! standard timevalue is seconds since ...
  integer, parameter  ::  since_time6(6) = (/1900,01,01,00,00,00/)


  !--- type --------------------------------------
  
  type TMeteoFile_tm5_nc
    ! input/output ?
    character(len=1)       ::  io
    !
    ! field collection
    !
    character(len=256)     ::  fname
    character(len=256)     ::  paramkeys      ! -aa-bb-cc-
    integer                ::  nparam
    character(len=64)      ::  tres
    type(TDate)            ::  trange(2)
    logical                ::  is_aver
    integer                ::  rnk
    !
    ! file
    !
    integer                ::  hid
    integer                ::  dimid_nv
    integer                ::  dimid_lon, varid_lon, varid_lon_bounds
    integer                ::  dimid_lat, varid_lat, varid_lat_bounds
    integer                ::  dimid_lonb, varid_lonb
    integer                ::  dimid_latb, varid_latb
    integer                ::  dimid_lev, varid_lev, varid_ap, varid_b, varid_ap_bounds, varid_b_bounds
    integer                ::  dimid_lev_u, varid_lev_u
    integer                ::  dimid_lev_v, varid_lev_v
    integer                ::  dimid_time, varid_time, varid_time_bounds, varid_reftime
    integer                ::  dimid_timeval, varid_timevalues, varid_timevalues_bounds, varid_reftimevalues
    integer                ::  varid_cell_area
    integer                ::  varid_ps, varid_ps_u, varid_ps_v
    !integer                ::  dimid_tm5_lm, dimid_tm5_lmb, varid_tm5_at, varid_tm5_bt
    integer, allocatable              ::  varid_param(:)
    character(len=16), allocatable    ::  varname_param(:)      ! (/'aa','bb','cc'/)
    character(len=256), allocatable   ::  cfname_param(:)       ! from CF table
    character(len=64), allocatable    ::  cfunit_param(:)       ! from CF table, following UDUnits
    integer, allocatable              ::  itrec_param(:)
    !
    ! output
    !
    logical                ::  output_initialised
    integer                ::  output_nrec
    integer                ::  output_ntrec
    !
    ! adhoc ...
    integer                ::  fixyear
  end type TMeteoFile_tm5_nc
  
  
  ! --- 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
#ifdef with_parallel_io_meteo
     module procedure mf_ReadRecord_parallel_io
#else
     module procedure mf_ReadRecord_serial_io
#endif
  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_NC_Init( rcf, status )
  
    use GO, only : TrcFile
    use GO, only : GO_Timer_Def
    use TMM_CF, only : TMM_CF_Init

    ! --- in/out ---------------------------------
    
    type(TRcFile), intent(in)           ::  rcf
    integer, intent(out)                ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_MF_TM5_NC_Init'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------
    
    ! init cf table and udunits package:
    call TMM_CF_Init( rcf, status )
    IF_NOTOK_RETURN(status=1)
    
    ! define timers:
    call GO_Timer_Def( itim_readrecord, 'tmm tm5-nc readrecord', status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TMM_MF_TM5_NC_Init


  ! ***


  subroutine TMM_MF_TM5_NC_Done( status )
  
    use TMM_CF, only : TMM_CF_Done

    ! --- in/out ---------------------------------
    
    integer, intent(out)                ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/TMM_MF_TM5_NC_Done'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------
    
    ! done with standard names table and udunits package:
    call TMM_CF_Done( status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TMM_MF_TM5_NC_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_nc), 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
    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
    
    ! --- begin --------------------------------

    ! store i/o :
    mf%io = io
    
    ! default flags:
    mf%is_aver = .false.
    mf%rnk     = 2         ! most 2D fields

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


    !
    ! main file
    !

    ! * set mf_filekey (uvsp,t,etc) and parmeters:
    select case ( paramkey )
      case ( 'sp', 'tsp' )
        mf_filekey = paramkey
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%nparam    = 1
      case ( 'mfu', 'mfv' )
        mf_filekey = 'mfuv'
        mf%paramkeys = '-mfu-mfv-'
        mf%nparam    = 2
        mf%rnk       = 3
      case ( 'mfw' )
        mf_filekey = 'mfw'
        mf%paramkeys = '-mfw-'
        mf%nparam    = 1
        mf%rnk       = 3
      case ( 'T' )
        mf_filekey   = 't'
        mf%paramkeys = '-T-'
        mf%nparam    = 1
        mf%rnk       = 3
      case ( 'Q' )
        mf_filekey   = 'q'
        mf%paramkeys = '-Q-'
        mf%nparam    = 1
        mf%rnk       = 3
      case ( 'CLWC', 'CIWC', 'CC', 'CCO', 'CCU', &
             'clwc', 'ciwc', 'cc', 'cco', 'ccu'  )
        mf_filekey   = 'cld'
        mf%paramkeys = '-CLWC-CIWC-CC-CCO-CCU-'
        mf%nparam    = 5
        mf%rnk       = 3
      case ( 'eu', 'ed', 'du', 'dd' ) ! computed online:  'cloud_base', 'cloud_top', 'cloud_lfs' 
        mf_filekey   = 'sub'
        !mf%paramkeys = '-eu-ed-du-dd-cloud_base-cloud_top-cloud_lfs-'
        mf%paramkeys = '-eu-ed-du-dd-'
        mf%nparam    = 4
        mf%rnk       = 3
      case ( 'K' )
        mf_filekey   = 'diffus'
        mf%paramkeys = '-K-'
        mf%nparam    = 1
        mf%rnk       = 3
      ! o constant fields
      case ( 'oro', 'lsm' )
        mf%tres      = 'constant'
        mf_filekey   = trim(paramkey)
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%nparam    = 1
      ! o monthly fields:
      case ( 'srols' )
        mf%tres      = 'month'
        mf_filekey   = trim(paramkey)
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%nparam    = 1
      ! o vegetation fields
      case ( 'cvl', 'cvh', &
             'tv01', 'tv02', 'tv03', 'tv04', 'tv05', &
             'tv06', 'tv07',         'tv09', 'tv10', &
             'tv11',         'tv13',                 &
             'tv16', 'tv17', 'tv18', 'tv19'         )
        mf_filekey = 'veg'
        mf%paramkeys = '-'//&
             'cvl-cvh-'//&
             'tv01-tv02-tv03-tv04-tv05-'//&
             'tv06-tv07-tv09-tv10-'//&
             'tv11-tv13-'//&
             'tv16-tv17-tv18-tv19-'
        mf%nparam    = 17
      ! o each surface file in a seperate file:
      case ( 'sr', 'srmer', &
             'swvl1', &
             'albedo', 'lsrh', 'ci', 'fg10','g10m', 'u10m', 'v10m', 'sd', &
             'blh', &
             't2m', 'd2m', &
             'sstr', 'src', 'raero', 'ustar', &
             'sst', 'sps', 'skt', 'ch4fire' )
        mf_filekey = trim(paramkey)
        mf%paramkeys = '-'//trim(paramkey)//'-'
        mf%nparam    = 1
      case ( 'sshf', 'slhf', 'ewss', 'nsss', 'lsp', 'cp', 'sf', &
             'ssr', 'ssrd', 'str', 'strd' )
          mf_filekey   = trim(paramkey)
          mf%paramkeys = '-'//trim(paramkey)//'-'
          mf%nparam    = 1
          mf%is_aver = .true.
      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
    
    ! storage:
    allocate( mf%varid_param  (mf%nparam) )
    allocate( mf%varname_param(mf%nparam) )
    allocate( mf%cfname_param (mf%nparam) )
    allocate( mf%cfunit_param (mf%nparam) )
    allocate( mf%itrec_param  (mf%nparam) )

    ! 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,'(6a)') trim(dir), mf_psep, trim(mf_mdir), mf_nsep, trim(mf_filekey), trim(output_ext)
      case ( 'month' )
        ! filename without day and forecast key:
        write (mf%fname,'(5a,"_",i4.4,i2.2,a)') &
                trim(dir), mf_psep, trim(mf_mdir), mf_nsep, trim(mf_filekey), ccyy, mm, trim(output_ext)
      case default
        ! filename including date:
        write (mf%fname,'(5a,"_",i4.4,2i2.2,3a)') &
                trim(dir), mf_psep, trim(mf_mdir), mf_nsep, &
                trim(mf_filekey), ccyy, mm, dd, trim(mf_fckey), trim(mf%tres), trim(output_ext)
    end select

    ! in case of output, not initialised yet ...
    mf%output_initialised = .false.

    ! 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 MDF, only : MDF_Close

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

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

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/mf_Done'
    
    ! --- begin -------------------------------------
    
    ! file opend for output ?
    if ( mf%output_initialised ) then
      ! close file:
      call MDF_Close( mf%hid, status )
      if ( status /= 0 ) then
        write (gol,'("closing file : ",a)') trim(mf%fname); call goErr
        TRACEBACK; status=1; return
      endif
    end if
    
    ! clear:
    deallocate( mf%varid_param   )
    deallocate( mf%varname_param )
    deallocate( mf%cfname_param  )
    deallocate( mf%cfunit_param  )
    deallocate( mf%itrec_param   )
    
    ! 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_nc), 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
  !  o nrec    :  number of time records in completed 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+96, 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 15 18 21   :  nrec=8
              !  f4  [ 96,120)  : 00 03 06 09 12    18      :  nrec=6
              !  f5  [120,144)  : 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 < 96.0 ) then
                nrec = 8
              else if ( dhr < 120 ) 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
  ! ***
  ! ******************************************************************

  
  !
  ! initialiase grid info from sds
  !
  
  subroutine lli_Init_mf( lli, nuv, mf, status )
  
    use Grid, only : TllGridInfo,   Init
    use MDF,  only : MDF_Inq_DimID, MDF_Inquire_Dimension
    use MDF,  only : MDF_Inq_VarID, MDF_Get_Var

    ! --- in/out ----------------------------------
    
    type(TllGridInfo),       intent(inout) ::  lli
    character(len=1),        intent(in)    ::  nuv
    type(TMeteoFile_tm5_nc), intent(in)    ::  mf
    integer,                 intent(out)   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/lli_Init_mf'
    
    ! --- local -----------------------------------
    
    integer             ::  dimid, varid
    real, allocatable   ::  values(:)
    real                ::  lon_deg, dlon_deg
    integer             ::  nlon
    real                ::  lat_deg, dlat_deg
    integer             ::  nlat
    
    ! --- begin ------------------------------------
    
    ! number of longitudes:
    call MDF_Inq_DimID( mf%hid, 'lon', dimid, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Inquire_Dimension( mf%hid, dimid, status, length=nlon )
    IF_NOTOK_RETURN(status=1)
    ! lon axis:
    allocate( values(nlon) )
    ! read:
    call MDF_Inq_VarID( mf%hid, 'lon', varid, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Get_Var( mf%hid, varid, values, status )
    IF_NOTOK_RETURN(status=1)
    ! extract:
    lon_deg  = values(1)
    dlon_deg = values(2) - values(1)
    ! clear:
    deallocate( values )

    ! number of latgitudes:
    call MDF_Inq_DimID( mf%hid, 'lat', dimid, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Inquire_Dimension( mf%hid, dimid, status, length=nlat )
    IF_NOTOK_RETURN(status=1)
    ! lat axis:
    allocate( values(nlat) )
    ! read:
    call MDF_Inq_VarID( mf%hid, 'lat', varid, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Get_Var( mf%hid, varid, values, status )
    IF_NOTOK_RETURN(status=1)
    ! extract:
    lat_deg  = values(1)
    dlat_deg = values(2) - values(1)
    ! clear:
    deallocate( values )
    
    ! define grid:    
    call Init( lli, lon_deg, dlon_deg, nlon, &
                    lat_deg, dlat_deg, nlat, status )
    IF_NOTOK_RETURN(status=1)
    
    ! ok
    status = 0
    
  end subroutine lli_Init_mf
  
  
  ! ***
  

  !
  ! initialiase level info from sds
  !
  
  subroutine levi_Init_mf( levi, mf, status )
  
    use Grid, only : TLevelInfo, Init
    use MDF, only : MDF_Inq_DimID, MDF_Inquire_Dimension
    use MDF, only : MDF_Inq_VarID, MDF_Get_Var
    use MDF, only : MDF_Get_Att, MDF_GLOBAL

    ! --- in/out ----------------------------------
    
    type(TLevelInfo), intent(out)           ::  levi
    type(TMeteoFile_tm5_nc), intent(in)     ::  mf
    integer, intent(out)                    ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/levi_Init_mf'
    
    ! --- local -----------------------------------
    
    integer                ::  dimid, varid
    integer                ::  lev
    character(len=1)       ::  nw
    integer                ::  lm
    real, allocatable      ::  at(:), bt(:)
    real, allocatable      ::  ap_bounds(:,:), b_bounds(:,:)
    
    ! --- begin ------------------------------------
    
    ! 2D or 3D ?
    select case ( mf%rnk )
    
      case ( 2 )

        ! set dummy values ...
        call Init( levi, 1, (/0.0,0.0/), (/0.0,0.0/), status )
        IF_NOTOK_RETURN(status=1)

      case ( 3 )

        ! read number of levels from global attribute:
        call MDF_Inq_DimID( mf%hid, 'lev', dimid, status )
        IF_NOTOK_RETURN(status=1)
        call MDF_Inquire_Dimension( mf%hid, dimid, status, length=lev )
        IF_NOTOK_RETURN(status=1)
        
        ! layers or levels ?
        call MDF_Get_Att( mf%hid, MDF_GLOBAL, 'nw', nw, status )
        IF_NOTOK_RETURN(status=1)
        
        ! extract:
        select case ( nw )
          !~ layers
          case ( 'n' )!, '*' )
            ! layers is lev ...
            lm = lev
            ! storage:
            allocate( at(lm+1), bt(lm+1) )
            allocate( ap_bounds(2,lev), b_bounds(2,lev) )
            ! extract hybride coeff
            call MDF_Inq_VarID( mf%hid, 'ap_bounds', varid, status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Get_Var( mf%hid, varid, ap_bounds, status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Inq_VarID( mf%hid, 'b_bounds', varid, status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Get_Var( mf%hid, varid, b_bounds, status )
            IF_NOTOK_RETURN(status=1)
            ! extract values:
            at = (/ ap_bounds(1,1), ap_bounds(2,:) /)
            bt = (/  b_bounds(1,1),  b_bounds(2,:) /)
            ! clear:
            deallocate( ap_bounds, b_bounds )
          !~ half levels
          case ( 'w' )
            ! layers is one less than lev ...
            lm = lev-1
            ! storage:
            allocate( at(lm+1), bt(lm+1) )
            ! extract hybride coeff
            call MDF_Inq_VarID( mf%hid, 'ap', varid, status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Get_Var( mf%hid, varid, at, status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Inq_VarID( mf%hid, 'b', varid, status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Get_Var( mf%hid, varid, bt, status )
            IF_NOTOK_RETURN(status=1)
          !~ unknown ...
          case default
            write (gol,'("found unsupported value of nw attribute : ",a)') trim(nw); call goErr
            write (gol,'("  file : ",a)') trim(mf%fname); call goErr
            TRACEBACK; status=1; return
        end select

        ! fill ...
        call Init( levi, lm, at, bt, status )
        IF_NOTOK_RETURN(status=1)

        ! clear:    
        deallocate( at, bt )

      case default
      
        write (gol,'("unsupported rank : ",i1)') mf%rnk; call goErr
        TRACEBACK; status=1; return
        
    end select
    
    ! ok
    status = 0
    
  end subroutine levi_Init_mf
    

  ! ***

  subroutine GetTimeRecordNumber( mf, t1, t2, irec, status )

    use GO, only : TDate, Get
    use MDF, only : MDF_Inq_DimID, MDF_Inquire_Dimension
    use MDF, only : MDF_Inq_VarID, MDF_Inquire_Variable, MDF_Get_Var

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), intent(inout)  ::  mf
    type(TDate), intent(in)                 ::  t1, t2
    integer, intent(out)                    ::  irec
    integer, intent(out)                    ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/GetTimeRecordNumber'
    
    ! --- local -------------------------------
    
    integer               ::  timevalues1(6), timevalues2(6)
    integer               ::  dimid, varid
    integer               ::  ntime, itime
    integer, allocatable  ::  timevalues_bounds(:,:,:)
    
    ! --- begin ---------------------------------
  
    ! target times:
    call Get( t1, time6=timevalues1 )
    call Get( t2, time6=timevalues2 )

    ! number of time records:
    call MDF_Inq_DimID( mf%hid, 'time', dimid, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Inquire_Dimension( mf%hid, dimid, status, length=ntime )
    IF_NOTOK_RETURN(status=1)

    ! storage:
    allocate( timevalues_bounds(2,6,ntime) )

    ! read time boundaries:
    call MDF_Inq_VarID( mf%hid, 'timevalues_bounds', varid, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Get_Var( mf%hid, varid, timevalues_bounds, status )
    IF_NOTOK_RETURN(status=1)

    ! search ...
    irec = -1
    do itime = 1, ntime
      ! compare:
      if ( all(timevalues_bounds(1,:,itime) == timevalues1) .and. &
           all(timevalues_bounds(2,:,itime) == timevalues2) ) then
        irec = itime
        exit
      end if
    end do

    ! check ...
    if ( irec < 0 ) then
      write (gol,'("could not find time record for:")'); call goErr
      write (gol,'("  t1, t2  : ",i4,2("-",i2.2)," ",i2.2,2(":",i2.2)," , ",'// &
                               & 'i4,2("-",i2.2)," ",i2.2,2(":",i2.2))') timevalues1, timevalues2; call goErr
      write (gol,'("in time bounds:")'); call goErr
      do itime = 1, ntime
        write (gol,'("  ",i6,"  : ",i4,2("-",i2.2)," ",i2.2,2(":",i2.2)," , ",'// &
                                 & 'i4,2("-",i2.2)," ",i2.2,2(":",i2.2))') &
                       itime, timevalues_bounds(1,:,itime), timevalues_bounds(2,:,itime); call goErr
      end do
      write (gol,'("in file:")'); call goErr
      write (gol,'("  ",a)') trim(mf%fname); call goErr
      TRACEBACK; status=1; return
    end if
    
    ! clear:
    deallocate( timevalues_bounds )

    ! ok
    status = 0

  end subroutine GetTimeRecordNumber
  
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    MF_READRECORD_PARALLEL_IO
  !
  ! !DESCRIPTION: Reads a met field in parallel. If the met field is 3D, also
  !                reads the surface pressure.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE MF_READRECORD_PARALLEL_IO( mf, paramkey, unit, t1, t2, nuv, nw,  &
                                            gridtype, levi, lli, ll, sp_ll, status )
    !
    ! !USES:
    !
    use PArray,   only : pa_Done, pa_SetShape
    use GO,       only : TDate
    use GO,       only : GO_Timer_Start, GO_Timer_End
    use Grid,     only : TllGridInfo, TLevelInfo, Done, operator(==)
    
    use MDF,      only : MDF_Open,      MDF_READ,              MDF_Close
    use MDF,      only : MDF_Inq_DimID, MDF_Inquire_Dimension, MDF_Var_Par_Access
    use MDF,      only : MDF_Inq_VarID, MDF_Inquire_Variable,  MDF_Get_Var
    use MDF,      only : MDF_Get_Att,   MDF_COLLECTIVE

    use partools, only : localComm,  MPI_INFO_NULL
    use dims,     only : nregions_all
    use TMM_CF,   only : TMM_CF_Convert_Units
    
    use tm5_distgrid, only : dgrid, dist_grid, Get_DistGrid, myid, Init_DistGrid, Done_DistGrid
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TMeteoFile_tm5_nc), intent(inout) :: mf
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*),        intent(in)    :: paramkey ! name of wanted variable
    character(len=*),        intent(in)    :: unit     ! wanted unit
    type(TDate),             intent(in)    :: t1, t2   ! time interval where to look for a time record
    character(len=1),        intent(in)    :: nuv      ! horizontal staggering
    character(len=1),        intent(in)    :: nw       ! vertical staggering
    !
    ! !OUTPUT PARAMETERS:
    !
    character(len=2),        intent(out)   :: gridtype   ! grid type of the  met field
    type(TLevelInfo),        intent(out)   :: levi       ! level info for METEO grid
    type(TllGridInfo),       intent(inout) :: lli        ! horizontal grid info of METEO grid
    real, pointer                          :: ll(:,:,:)  ! read met field 
    real, pointer                          :: sp_ll(:,:) ! read surface pressure
    integer,                 intent(out)   :: status     ! error code
    !
    ! !REVISION HISTORY: 
    !   17 Oct 2013 - Ph. Le Sager - v0
    !
    ! !REMARKS:
    ! - (1) Cannot use the following because it creates a circular dependency:
    !    use meteodata,    only : global_lli
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/mf_ReadRecord_parallel_io'

    integer           :: varid, klm
    integer           :: irec
    integer           :: shp(7), ndims, access_mode, sz(3)
    character(len=64) :: cfunits
    real              :: ufac
    integer           :: i1, i2, j1, j2      ! bounds of the tile w/r/t the METEO grid!
    type(TllGridInfo) :: world_lli
    type(dist_grid)   :: DistGrid
    
    ! --- 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

    ! access mode & open
    access_mode = MDF_COLLECTIVE

    call MDF_Open( trim(mf%fname), MDF_NETCDF4, MDF_READ, mf%hid, status, mpi_comm=localComm, mpi_info=MPI_INFO_NULL)
    IF_NOTOK_RETURN(status=1)
    

    ! *** VARIABLE ID

    call MDF_Inq_VarID( mf%hid, trim(paramkey), varid, status )
    IF_NOTOK_RETURN(status=1)

    call MDF_Var_Par_Access( mf%hid, varid, access_mode, status )
    IF_NOTOK_RETURN(status=1)


    ! *** GRID DEFINITION

    gridtype = 'll'

    ! Global grid info: Find from file's sds and init
    call lli_Init_mf( lli, nuv, mf, status )
    IF_NOTOK_RETURN(status=1)

    ! setup level definition
    call levi_Init_mf( levi, mf, status )
    IF_NOTOK_RETURN(status=1)

    ! get global dimensions (here we are interested by the third dimension only)
    call MDF_Inquire_Variable( mf%hid, varid, status, ndims=ndims )
    IF_NOTOK_RETURN(status=1)
    
    call MDF_Inquire_Variable( mf%hid, varid, status, shp=shp(1:ndims) )
    IF_NOTOK_RETURN(status=1)

    ! Get distributed grid information.
    ! Assuming that the meteo grid is also one of the 1..nregions_all, find which one:
    do klm=1, nregions_all
       call Get_DistGrid( dgrid(klm), global_lli=world_lli )

       if (lli == world_lli) exit
    enddo
    call done(world_lli, status)     ! garbage clean
    IF_NOTOK_RETURN(status=1)
    
    ! Get local bounds and local lli, from either an existing dgrid, or from a
    ! newly created one if no match
    if ( klm == (nregions_all+1) ) then  ! no match

       call Init_DistGrid( DistGrid, lli, myid, 0, status)
       IF_NOTOK_RETURN(status=1)
       
       call Get_DistGrid( DistGrid, I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, lli=lli )

       call Done_DistGrid(DistGrid, status)
       IF_NOTOK_RETURN(status=1)
       
    else
       
       call Get_DistGrid( dgrid(klm), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2, lli=lli )
    end if

    
    ! *** WORK ARRAY

    ! shape to read - Note overlap between tiles
    select case ( nuv )
    case ( 'n' )
       sz=(/i2-i1+1,j2-j1+1,shp(3)/)
    case ( 'u' )
       sz=(/i2-i1+2,j2-j1+1,shp(3)/)
    case ( 'v' )
       sz=(/i2-i1+1,j2-j1+2,shp(3)/)
    end select
       

    ! *** EXTRACT DATA

    if ( trim(mf%tres) == 'constant' ) then      ! constant field

       select case ( ndims )

       case ( 2 )
          allocate( ll(sz(1),sz(2),1) )

          call MDF_Get_Var( mf%hid, varid, ll(:,:,1), status, start=(/i1,j1/), count=sz(1:2) )
          IF_NOTOK_RETURN(status=1)

       case ( 3 )
          call pa_SetShape( ll, sz )

          call MDF_Get_Var( mf%hid, varid, ll, status, start=(/i1,j1,1/), count=sz )
          IF_NOTOK_RETURN(status=1)

       case default
          write (gol,'("unsupported data rank:",i6)') ndims; call goErr
          TRACEBACK; status=1; return
       end select

    else

       ! which time record ?
       call GetTimeRecordNumber( mf, t1, t2, irec, status )
       IF_NOTOK_RETURN(status=1)

       select case ( ndims )

       case ( 2+1 )
          allocate( ll(sz(1),sz(2),1) )

          call MDF_Get_Var( mf%hid, varid, ll(:,:,1), status, start=(/i1,j1,irec/), count=(/sz(1),sz(2),1/) )
          IF_NOTOK_RETURN(status=1)

       case ( 3+1 )
          call pa_SetShape( ll, sz )

          call MDF_Get_Var( mf%hid, varid, ll, status, start=(/i1,j1,1,irec/), count=(/sz(1),sz(2),sz(3),1/) )
          IF_NOTOK_RETURN(status=1)

       case default
          write (gol,'("unsupported data rank:",i6)') ndims; call goErr
          TRACEBACK; status=1; return
       end select

    end if

    
    ! *** UNIT CONVERSION

    ! get unit from field in file:
    call MDF_Get_Att( mf%hid, varid, 'units', cfunits, status )
    IF_NOTOK_RETURN(status=1)

!PLS    ! conversion factor:
!PLS    call TMM_CF_Convert_Units( trim(cfunits), trim(unit), ufac, status )
!PLS    IF_NOTOK_RETURN(status=1)
!PLS    ! apply ?
!PLS    if ( ufac /= 1.0 ) then
!PLS       ! convert:
!PLS       ll = ll * ufac
!PLS       !! info ...
!PLS       !write (gol,'("      convert `",a,"` from `",a,"` to `",a,"` with factor ",f8.2," ; new range ",2f8.2)') &
!PLS       !               trim(paramkey), trim(cfunits), trim(unit), ufac, minval(ll), maxval(ll); call goPr
!PLS    end if

    
    ! *** SURFACE PRESSURE

    if ( mf%rnk == 3 ) then  ! required?
       ! search variable:
       select case ( nuv )
       case ( 'n' )
          call MDF_Inq_VarID( mf%hid, 'ps', varid, status )
          IF_NOTOK_RETURN(status=1)
       case ( 'u' )
          call MDF_Inq_VarID( mf%hid, 'ps_u', varid, status )
          IF_NOTOK_RETURN(status=1)
       case ( 'v' )
          call MDF_Inq_VarID( mf%hid, 'ps_v', varid, status )
          IF_NOTOK_RETURN(status=1)
       case default
          write (gol,'("unsupported nuv :",a)') nuv; call goErr
          TRACEBACK; status=1; return
       end select

       call MDF_Var_Par_Access( mf%hid, varid, access_mode, status )
       IF_NOTOK_RETURN(status=1)

       call pa_SetShape( sp_ll, sz(1:2) )

       call MDF_Get_Var( mf%hid, varid, sp_ll, status, start=(/i1,j1,irec/), count=(/sz(1),sz(2),1/) )
       IF_NOTOK_RETURN(status=1)
    else
       ! for safety ...
       call pa_Done( sp_ll )
    end if

    
    ! *** DONE

    call MDF_Close( mf%hid, status )
    IF_NOTOK_RETURN(status=1)

    call GO_Timer_End( itim_readrecord, status )
    IF_NOTOK_RETURN(status=1)

    status = 0

  END SUBROUTINE MF_READRECORD_PARALLEL_IO
  !EOC

  

  subroutine mf_ReadRecord_serial_io( mf, paramkey, unit, t1, t2, nuv, nw, &
                                          gridtype, levi, &
                                          lli, ll, sp_ll, &
                                          status )

    use PArray, only : pa_Done, pa_SetShape
    use GO, only : TDate
    use GO, only : GO_Timer_Start, GO_Timer_End
    use Grid, only : TllGridInfo, TLevelInfo
    use MDF, only : MDF_Open, MDF_Close
    use MDF, only : MDF_READ
    use MDF, only : MDF_Inq_DimID, MDF_Inquire_Dimension
    use MDF, only : MDF_Inq_VarID, MDF_Inquire_Variable, MDF_Get_Var
    use MDF, only : MDF_Get_Att

    use TMM_CF, only : TMM_CF_Convert_Units

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), intent(inout) ::  mf
    character(len=*),        intent(in)    ::  paramkey
    character(len=*),        intent(in)    ::  unit
    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 -------------------------------
    
    integer               ::  varid
    integer               ::  irec
    integer               ::  shp(7), ndims
    character(len=64)     ::  cfunits
    real                  ::  ufac
    
    ! --- 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:
    call MDF_Open( trim(mf%fname), output_type, MDF_READ, mf%hid, status )
    IF_NOTOK_RETURN(status=1)
    
    ! *** variable id
    
    ! search variable:
    call MDF_Inq_VarID( mf%hid, trim(paramkey), varid, status )
    IF_NOTOK_RETURN(status=1)

    ! *** grid definition
    
    ! always regular lat/lon grid ..    
    gridtype = 'll'
    
    ! setup grid definition:
    call lli_Init_mf( lli, nuv, mf, status )
    IF_NOTOK_RETURN(status=1)
    
    ! setup level definition:
    call levi_Init_mf( levi, mf, status )
    IF_NOTOK_RETURN(status=1)
    
    ! get dimensions:
    call MDF_Inquire_Variable( mf%hid, varid, status, ndims=ndims )
    IF_NOTOK_RETURN(status=1)
    call MDF_Inquire_Variable( mf%hid, varid, status, shp=shp(1:ndims) )
    IF_NOTOK_RETURN(status=1)
    
    ! *** data
    
    ! constant field ?
    if ( trim(mf%tres) == 'constant' ) then

      ! extract data array:
      select case ( ndims )
        !* 2D fields
        case ( 2 )
          ! storage:
          shp(3) = 1
          call pa_SetShape( ll, shp(1:3) )
          ! complete record:
          call MDF_Get_Var( mf%hid, varid, ll(:,:,1), status )
          IF_NOTOK_RETURN(status=1)
        !* 3d fields
        case ( 3 )
          ! storage:
          call pa_SetShape( ll, shp(1:3) )
          ! complete record:
          call MDF_Get_Var( mf%hid, varid, ll, status )
          IF_NOTOK_RETURN(status=1)
        !* unknown
        case default
          write (gol,'("unsupported data rank:",i6)') ndims; call goErr
          TRACEBACK; status=1; return
      end select
    
    else

      ! which time record ?
      call GetTimeRecordNumber( mf, t1, t2, irec, status )
      IF_NOTOK_RETURN(status=1)

      ! extract data array:
      select case ( ndims )
        !* 2D temporal fields
        case ( 2+1 )
          ! storage:
          shp(3) = 1
          call pa_SetShape( ll, shp(1:3) )
          ! complete record:
          call MDF_Get_Var( mf%hid, varid, ll(:,:,1), status, &
                                start=(/1,1,irec/), count=(/shp(1),shp(2),1/) )
          IF_NOTOK_RETURN(status=1)
        !* 3d temporal fields
        case ( 3+1 )
          ! storage:
          call pa_SetShape( ll, shp(1:3) )
          ! complete record:
          call MDF_Get_Var( mf%hid, varid, ll, status, &
                                start=(/1,1,1,irec/), count=(/shp(1),shp(2),shp(3),1/) )
          IF_NOTOK_RETURN(status=1)
        !* unknown
        case default
          write (gol,'("unsupported data rank:",i6)') ndims; call goErr
          TRACEBACK; status=1; return
      end select
      
    end if
    
    ! *** unit conversion

    ! get unit from field in file:
    call MDF_Get_Att( mf%hid, varid, 'units', cfunits, status )
    IF_NOTOK_RETURN(status=1)
!PLS   ! conversion factor:
!PLS   call TMM_CF_Convert_Units( trim(cfunits), trim(unit), ufac, status )
!PLS   IF_NOTOK_RETURN(status=1)
!PLS   ! apply ?
!PLS   if ( ufac /= 1.0 ) then
!PLS     ! convert:
!PLS     ll = ll * ufac
!PLS     !! info ...
!PLS     !write (gol,'("      convert `",a,"` from `",a,"` to `",a,"` with factor ",f8.2," ; new range ",2f8.2)') &
!PLS     !               trim(paramkey), trim(cfunits), trim(unit), ufac, minval(ll), maxval(ll); call goPr
!PLS   end if
    
    ! *** surface pressure

    ! surface pressure required ?
    if ( mf%rnk == 3 ) then
      ! search variable:
      select case ( nuv )
        case ( 'n' )
          call MDF_Inq_VarID( mf%hid, 'ps', varid, status )
          IF_NOTOK_RETURN(status=1)
        case ( 'u' )
          call MDF_Inq_VarID( mf%hid, 'ps_u', varid, status )
          IF_NOTOK_RETURN(status=1)
        case ( 'v' )
          call MDF_Inq_VarID( mf%hid, 'ps_v', varid, status )
          IF_NOTOK_RETURN(status=1)
        case default
          write (gol,'("unsupported nuv :",a)') nuv; call goErr
          TRACEBACK; status=1; return
      end select
      ! storage:
      call pa_SetShape( sp_ll, shp(1:2) )
      ! complete record:
      call MDF_Get_Var( mf%hid, varid, sp_ll, status, &
                            start=(/1,1,irec/), count=(/shp(1),shp(2),1/) )
      IF_NOTOK_RETURN(status=1)
    else
      ! for safety ...
      call pa_Done( sp_ll )
    end if
    
    ! ***
 
    ! close
    call MDF_Close( mf%hid, 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_serial_io
  

  ! ******************************************************************
  ! ***
  ! *** output
  ! ***
  ! ******************************************************************
  
  
  subroutine CF_Put_Standard_Atts( hid, varid, cf_standard_name, cf_units, status )
  
    use MDF, only : MDF_Put_Att
    
    ! --- in/out ---------------------------------
    
    integer, intent(in)                           ::  hid
    integer, intent(in)                           ::  varid
    character(len=*), intent(in)                  ::  cf_standard_name
    character(len=*), intent(in)                  ::  cf_units
    integer, intent(out)                          ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/CF_Put_Standard_Atts'
    
    ! --- local ----------------------------------
    
    ! --- begin ----------------------------------

    ! add standard name attribute:
    call MDF_Put_Att( hid, varid, 'standard_name', trim(cf_standard_name), status )
    IF_NOTOK_RETURN(status=1)
    
    ! add units attribute:
    call MDF_Put_Att( hid, varid, 'units', trim(cf_units), status )
    IF_NOTOK_RETURN(status=1)
    
    ! ok
    status = 0
    
  end subroutine CF_Put_Standard_Atts


  ! ***
  
  
  subroutine Define_File( mf, lli, nuv, status, levi, nw )!, nlev )

    use GO    , only : goReadFromLine
    use Binas , only : grav, ae
    use Grid  , only : TllGridInfo, TLevelInfo, AreaOper
    use MDF   , only : MDF_Put_Att, MDF_GLOBAL
    use MDF   , only : MDF_Def_Var, MDF_Put_Var, MDF_FLOAT, MDF_DOUBLE, MDF_INT
    use MDF   , only : MDF_Def_Dim, MDF_UNLIMITED
    use MDF   , only : MDF_EndDef
    use TMM_CF, only : TMM_CF_Standard_Units, TMM_CF_Convert_Name

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), intent(inout)   ::  mf
    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//'/Define_File'
    
    ! --- local ----------------------------------
    
    character(len=512)              ::  units
    character(len=512)              ::  cell_measure
    character(len=512)              ::  cell_methods
    character(len=512)              ::  coordinates
    real, allocatable               ::  pat(:,:)
    
    integer                         ::  l, n
    character(len=1024)             ::  line
    character(len=32)               ::  paramkey
    integer                         ::  iparam
    character(len=1024)             ::  cf_standard_name
    character(len=512)              ::  cf_units
    integer                         ::  varid
    
    ! --- begin ----------------------------------
    
    ! 
    ! CF standard global attribues
    !
    
    ! convention id:
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'Conventions', 'CF-1.4', status )
    IF_NOTOK_RETURN(status=1)
    
    ! descriptive title:
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'title', 'TM meteo file', status )
    IF_NOTOK_RETURN(status=1)
    
    ! who?
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'institution', 'TM community', status )
    IF_NOTOK_RETURN(status=1)
    
    ! from ?
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'source', 'TM produced meteo file', status )
    IF_NOTOK_RETURN(status=1)
    
    ! how ?
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'history', 'None', status )
    IF_NOTOK_RETURN(status=1)
    
    ! published material:
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'references', 'None', status )
    IF_NOTOK_RETURN(status=1)
    
    ! other:
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'comment', 'None', status )
    IF_NOTOK_RETURN(status=1)
    
    !
    ! TMM specific global attributes
    !
    
    ! file format
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'tmm_format', trim(output_format), status )
    IF_NOTOK_RETURN(status=1)
    ! grid type
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'tmm_gridtype', 'll', status )
    IF_NOTOK_RETURN(status=1)
    ! gravity constant:
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'grav', grav, status )
    IF_NOTOK_RETURN(status=1)
    ! earth radius:
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'ae', ae, status )
    IF_NOTOK_RETURN(status=1)
    
    ! save first and last lon/lat (center) for use with HIPHOP
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'lonmin', lli%lon_deg(1)       , status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'lonmax', lli%lon_deg(lli%nlon), status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'latmin', lli%lat_deg(1)       , status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'latmax', lli%lat_deg(lli%nlat), status )
    IF_NOTOK_RETURN(status=1)

    !
    ! axis
    !

    ! extra coordinate axes (with other names than the dimensions):
    coordinates = ''
    
    ! auxilary cell info:
    cell_measure = ''
    
    ! how data was formed:
    cell_methods = ''

    ! vertices:
    call MDF_Def_Dim( mf%hid, 'nv', 2, mf%dimid_nv, status )
    IF_NOTOK_RETURN(status=1)
    
    ! * longitudes
    
    ! longitude dimension:
    call MDF_Def_Dim( mf%hid, 'lon', lli%nlon, mf%dimid_lon, status )
    IF_NOTOK_RETURN(status=1)
    
    ! variable:
    call MDF_Def_Var( mf%hid, 'lon', MDF_DOUBLE, (/mf%dimid_lon/), mf%varid_lon, status )
    IF_NOTOK_RETURN(status=1)

    ! convert to CF name and unit:
    cf_standard_name = 'longitude'
    call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
    IF_NOTOK_RETURN(status=1)
    ! add standard_name and units attributes:
    call CF_Put_Standard_Atts( mf%hid, mf%varid_lon, trim(cf_standard_name), trim(cf_units), status )
    IF_NOTOK_RETURN(status=1)

    ! generic axis name:
    call MDF_Put_Att( mf%hid, mf%varid_lon, 'axis', 'X', status )
    IF_NOTOK_RETURN(status=1)
    
    ! add attribute with name of variable with boundaries:
    call MDF_Put_Att( mf%hid, mf%varid_lon, 'bounds', 'lon_bounds', status )
    IF_NOTOK_RETURN(status=1)
    ! add variable for boundaries:
    call MDF_Def_Var( mf%hid, 'lon_bounds', MDF_DOUBLE, (/mf%dimid_nv,mf%dimid_lon/), mf%varid_lon_bounds, status )
    IF_NOTOK_RETURN(status=1)
    
    ! extra:
    if ( nuv /= 'n' ) then
      ! dimension:
      call MDF_Def_Dim( mf%hid, 'lonb', lli%nlon+1, mf%dimid_lonb, status )
      IF_NOTOK_RETURN(status=1)
      ! variable:
      call MDF_Def_Var( mf%hid, 'lonb', MDF_DOUBLE, (/mf%dimid_lonb/), mf%varid_lonb, status )
      IF_NOTOK_RETURN(status=1)
      ! add standard_name and units attributes:
      call CF_Put_Standard_Atts( mf%hid, mf%varid_lonb, trim(cf_standard_name), trim(cf_units), status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! * latitudes

    ! latitude dimension:
    call MDF_Def_Dim( mf%hid, 'lat', lli%nlat, mf%dimid_lat, status )
    IF_NOTOK_RETURN(status=1)
    
    ! variable:
    call MDF_Def_Var( mf%hid, 'lat', MDF_DOUBLE, (/mf%dimid_lat/), mf%varid_lat, status )
    IF_NOTOK_RETURN(status=1)

    ! convert to CF name and unit:
    cf_standard_name = 'latitude'
    call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
    IF_NOTOK_RETURN(status=1)
    ! add standard_name and units attributes:
    call CF_Put_Standard_Atts( mf%hid, mf%varid_lat, trim(cf_standard_name), trim(cf_units), status )
    IF_NOTOK_RETURN(status=1)

    ! generic axis name:
    call MDF_Put_Att( mf%hid, mf%varid_lat, 'axis', 'Y', status )
    IF_NOTOK_RETURN(status=1)
    
    ! add attribute with name of variable with boundaries:
    call MDF_Put_Att( mf%hid, mf%varid_lat, 'bounds', 'lat_bounds', status )
    IF_NOTOK_RETURN(status=1)
    ! add variable for boundaries:
    call MDF_Def_Var( mf%hid, 'lat_bounds', MDF_DOUBLE, (/mf%dimid_nv,mf%dimid_lat/), mf%varid_lat_bounds, status )
    IF_NOTOK_RETURN(status=1)
    
    ! extra:
    if ( nuv /= 'n' ) then
      ! dimension:
      call MDF_Def_Dim( mf%hid, 'latb', lli%nlat+1, mf%dimid_latb, status )
      IF_NOTOK_RETURN(status=1)
      ! variable:
      call MDF_Def_Var( mf%hid, 'latb', MDF_DOUBLE, (/mf%dimid_latb/), mf%varid_latb, status )
      IF_NOTOK_RETURN(status=1)
      ! add standard_name and units attributes:
      call CF_Put_Standard_Atts( mf%hid, mf%varid_latb, trim(cf_standard_name), trim(cf_units), status )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! * area

    ! cell variable ?
    if ( nuv == 'n' ) then
      ! also provide the area:
      cell_measure = trim(cell_measure)//' area: cell_area'
      ! cell area:
      call MDF_Def_Var( mf%hid, 'cell_area', MDF_FLOAT, (/mf%dimid_lon,mf%dimid_lat/), mf%varid_cell_area, status )
      IF_NOTOK_RETURN(status=1)
      ! convert to CF name and unit:
      cf_standard_name = 'cell_area'
      call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
      IF_NOTOK_RETURN(status=1)
      ! add standard_name and units attributes:
      call CF_Put_Standard_Atts( mf%hid, mf%varid_cell_area, trim(cf_standard_name), trim(cf_units), status )
      IF_NOTOK_RETURN(status=1)
      ! add description:
      call MDF_Put_Att( mf%hid, mf%varid_cell_area, 'long_name', 'area of grid cell', status )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! * time
    
    if ( trim(mf%tres) /= 'constant' ) then

      ! time dimension:
      call MDF_Def_Dim( mf%hid, 'time', MDF_UNLIMITED, mf%dimid_time, status )
      IF_NOTOK_RETURN(status=1)

      ! standard units is single value only ...
      write (units,'("seconds since ",i4.4,2("-",i2.2)," ",i2.2,2(":",i2.2))') since_time6

      ! variable:
      call MDF_Def_Var( mf%hid, 'time', MDF_DOUBLE, (/mf%dimid_time/), mf%varid_time, status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_time, 'standard_name', 'time', status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_time, 'units', trim(units), status )
      IF_NOTOK_RETURN(status=1)

      ! name of variable with interval bounds:
      call MDF_Put_Att( mf%hid, mf%varid_time, 'bounds', 'time_bounds', status )
      IF_NOTOK_RETURN(status=1)

      ! variable:
      call MDF_Def_Var( mf%hid, 'time_bounds', MDF_DOUBLE, (/mf%dimid_nv,mf%dimid_time/), mf%varid_time_bounds, status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_time, 'units', trim(units), status )
      IF_NOTOK_RETURN(status=1)

      ! time averages ?
      if ( mf%is_aver ) then
        ! value is the mean over a time interval:
        cell_methods = trim(cell_methods)//' time: mean'
      else
        ! instant fields:
        cell_methods = trim(cell_methods)//' time: point'
      end if

      ! variable:
      call MDF_Def_Var( mf%hid, 'reftime', MDF_DOUBLE, (/mf%dimid_time/), mf%varid_reftime, status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_reftime, 'units', trim(units), status )
      IF_NOTOK_RETURN(status=1)
      ! description:
      call MDF_Put_Att( mf%hid, mf%varid_reftime, 'long_name', 'reference time', status )
      IF_NOTOK_RETURN(status=1)
      ! extra time coordinates:
      coordinates = trim(coordinates)//' reftime'
      
    end if
    
    ! * time values
    
    if ( trim(mf%tres) /= 'constant' ) then
    
      ! time values dimension:
      call MDF_Def_Dim( mf%hid, 'timeval', 6, mf%dimid_timeval, status )
      IF_NOTOK_RETURN(status=1)

      ! variable:
      call MDF_Def_Var( mf%hid, 'timevalues', MDF_INT, (/mf%dimid_timeval,mf%dimid_time/), mf%varid_timevalues, status )
      IF_NOTOK_RETURN(status=1)
      ! description:
      call MDF_Put_Att( mf%hid, mf%varid_timevalues, 'long_name', 'year month day hour minute second', status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_timevalues, 'units', '1', status )
      IF_NOTOK_RETURN(status=1)

      ! name of variable with interval bounds:
      call MDF_Put_Att( mf%hid, mf%varid_timevalues, 'bounds', 'timevalues_bounds', status )
      IF_NOTOK_RETURN(status=1)

      ! variable:
      call MDF_Def_Var( mf%hid, 'timevalues_bounds', MDF_INT, &
                           (/mf%dimid_nv,mf%dimid_timeval,mf%dimid_time/), &
                           mf%varid_timevalues_bounds, status )
      IF_NOTOK_RETURN(status=1)
      ! description:
      call MDF_Put_Att( mf%hid, mf%varid_timevalues_bounds, 'long_name', 'year month day hour minute second', status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_timevalues_bounds, 'units', '1', status )
      IF_NOTOK_RETURN(status=1)

      ! variable:
      call MDF_Def_Var( mf%hid, 'reftimevalues', MDF_INT, (/mf%dimid_timeval,mf%dimid_time/), mf%varid_reftimevalues, status )
      IF_NOTOK_RETURN(status=1)
      ! description:
      call MDF_Put_Att( mf%hid, mf%varid_reftimevalues, 'long_name', 'reference time values', status )
      IF_NOTOK_RETURN(status=1)
      ! add units attribute:
      call MDF_Put_Att( mf%hid, mf%varid_reftimevalues, 'units', '1', status )
      IF_NOTOK_RETURN(status=1)
      
    end if
    
    ! * levels
    
    ! levels ?
    if ( present(levi) ) then
    
      ! check ...
      if ( .not. present(nw) ) then
        write (gol,'("optional argument levi without nw ...")'); call goErr
        TRACEBACK; status=1; return
      end if
      
      !! tm5 specials ...
      !call MDF_Def_Dim( mf%hid, 'tm5_lm', levi%nlev, mf%dimid_tm5_lm, status )
      !IF_NOTOK_RETURN(status=1)
      !call MDF_Def_Dim( mf%hid, 'tm5_lmb', levi%nlev+1, mf%dimid_tm5_lmb, status )
      !IF_NOTOK_RETURN(status=1)
      !call MDF_Def_Var( mf%hid, 'tm5_at', MDF_DOUBLE, (/mf%dimid_tm5_lmb/), mf%varid_tm5_at, status )
      !IF_NOTOK_RETURN(status=1)
      !call MDF_Def_Var( mf%hid, 'tm5_bt', MDF_DOUBLE, (/mf%dimid_tm5_lmb/), mf%varid_tm5_bt, status )
      !IF_NOTOK_RETURN(status=1)
      
      ! save nw to facilitate reading:
      call MDF_Put_Att( mf%hid, MDF_GLOBAL, 'nw', trim(nw), status )
      IF_NOTOK_RETURN(status=1)

      ! where defined ?
      select case ( nw )
      
        ! layer mid, or selection
        case ( 'n' )!, '*' )
    
          ! hybride layers:
          !if ( present(nlev) ) then
          !  ! lowest layers only:
          !  call MDF_Def_Dim( mf%hid, 'lev', nlev, mf%dimid_lev, status )
          !  IF_NOTOK_RETURN(status=1)
          !else
            ! full dimension:
            call MDF_Def_Dim( mf%hid, 'lev', levi%nlev, mf%dimid_lev, status )
            IF_NOTOK_RETURN(status=1)
          !end if

          ! variable:
          call MDF_Def_Var( mf%hid, 'lev', MDF_DOUBLE, (/mf%dimid_lev/), mf%varid_lev, status )
          IF_NOTOK_RETURN(status=1)
          ! convert to CF name and unit:
          cf_standard_name = 'atmosphere_hybrid_sigma_pressure_coordinate'
          call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
          IF_NOTOK_RETURN(status=1)
          ! add standard_name and units attributes:
          call CF_Put_Standard_Atts( mf%hid, mf%varid_lev, trim(cf_standard_name), trim(cf_units), status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'long_name', 'pressure at layer midpoints', status )
          IF_NOTOK_RETURN(status=1)
          ! direction of increasing pressure:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'positive', 'down', status )
          IF_NOTOK_RETURN(status=1)
          ! how to compute pressure:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'forumula', 'p(n,k,j,i) = ap(k) + b(k)*ps(n,j,i)', status )
          IF_NOTOK_RETURN(status=1)
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'forumula_terms', 'ap: ap b: b ps: ps', status )
          IF_NOTOK_RETURN(status=1)
          ! generic axis name:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'axis', 'Z', status )
          IF_NOTOK_RETURN(status=1)
          
          ! extra
          if ( nuv /= 'n' ) then
            ! hybride layers:
            call MDF_Def_Dim( mf%hid, 'lev_u', levi%nlev, mf%dimid_lev_u, status )
            IF_NOTOK_RETURN(status=1)
            ! variable:
            call MDF_Def_Var( mf%hid, 'lev_u', MDF_DOUBLE, (/mf%dimid_lev_u/), mf%varid_lev_u, status )
            IF_NOTOK_RETURN(status=1)
            ! convert to CF name and unit:
            cf_standard_name = 'atmosphere_hybrid_sigma_pressure_coordinate'
            call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
            IF_NOTOK_RETURN(status=1)
            ! add standard_name and units attributes:
            call CF_Put_Standard_Atts( mf%hid, mf%varid_lev_u, trim(cf_standard_name), trim(cf_units), status )
            IF_NOTOK_RETURN(status=1)
            ! description:
            call MDF_Put_Att( mf%hid, mf%varid_lev_u, 'long_name', 'pressure at layer midpoints', status )
            IF_NOTOK_RETURN(status=1)
            ! direction of increasing pressure:
            call MDF_Put_Att( mf%hid, mf%varid_lev_u, 'positive', 'down', status )
            IF_NOTOK_RETURN(status=1)
            ! how to compute pressure:
            call MDF_Put_Att( mf%hid, mf%varid_lev_u, 'forumula', 'p(n,k,j,i) = ap(k) + b(k)*ps(n,j,i)', status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Put_Att( mf%hid, mf%varid_lev_u, 'forumula_terms', 'ap: ap b: b ps: ps_u', status )
            IF_NOTOK_RETURN(status=1)

            ! hybride layers:
            call MDF_Def_Dim( mf%hid, 'lev_v', levi%nlev, mf%dimid_lev_v, status )
            IF_NOTOK_RETURN(status=1)
            ! variable:
            call MDF_Def_Var( mf%hid, 'lev_v', MDF_DOUBLE, (/mf%dimid_lev_v/), mf%varid_lev_v, status )
            IF_NOTOK_RETURN(status=1)
            ! convert to CF name and unit:
            cf_standard_name = 'atmosphere_hybrid_sigma_pressure_coordinate'
            call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
            IF_NOTOK_RETURN(status=1)
            ! add standard_name and units attributes:
            call CF_Put_Standard_Atts( mf%hid, mf%varid_lev_v, trim(cf_standard_name), trim(cf_units), status )
            IF_NOTOK_RETURN(status=1)
            ! description:
            call MDF_Put_Att( mf%hid, mf%varid_lev_v, 'long_name', 'pressure at layer midpoints', status )
            IF_NOTOK_RETURN(status=1)
            ! direction of increasing pressure:
            call MDF_Put_Att( mf%hid, mf%varid_lev_v, 'positive', 'down', status )
            IF_NOTOK_RETURN(status=1)
            ! how to compute pressure:
            call MDF_Put_Att( mf%hid, mf%varid_lev_v, 'forumula', 'p(n,k,j,i) = ap(k) + b(k)*ps(n,j,i)', status )
            IF_NOTOK_RETURN(status=1)
            call MDF_Put_Att( mf%hid, mf%varid_lev_v, 'forumula_terms', 'ap: ap b: b ps: ps_v', status )
            IF_NOTOK_RETURN(status=1)
          end if

          ! variable:
          call MDF_Def_Var( mf%hid, 'ap', MDF_DOUBLE, (/mf%dimid_lev/), mf%varid_ap, status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_ap, 'long_name', 'atmosphere hybrid sigma pressure coefficient ap&
               & at layer midpoints', status )
          IF_NOTOK_RETURN(status=1)
          ! units:
          call MDF_Put_Att( mf%hid, mf%varid_ap, 'units', 'Pa', status )
          IF_NOTOK_RETURN(status=1)
          ! name of variable with boundary values:
          call MDF_Put_Att( mf%hid, mf%varid_ap, 'bounds', 'ap_bounds', status )
          IF_NOTOK_RETURN(status=1)

          ! variable:
          call MDF_Def_Var( mf%hid, 'b', MDF_DOUBLE, (/mf%dimid_lev/), mf%varid_b, status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_b, 'long_name', 'atmosphere hybrid sigma pressure coefficient b&
               & at layer midpoints', status )
          IF_NOTOK_RETURN(status=1)
          ! units:
          call MDF_Put_Att( mf%hid, mf%varid_b, 'units', '1', status )
          IF_NOTOK_RETURN(status=1)
          ! name of variable with boundary values:
          call MDF_Put_Att( mf%hid, mf%varid_b, 'bounds', 'b_bounds', status )
          IF_NOTOK_RETURN(status=1)

          ! variable:
          call MDF_Def_Var( mf%hid, 'ap_bounds', MDF_DOUBLE, (/mf%dimid_nv,mf%dimid_lev/), mf%varid_ap_bounds, status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_ap_bounds, 'long_name', 'atmosphere hybrid sigma pressure coefficient&
               & ap at layer interfaces', status )
          IF_NOTOK_RETURN(status=1)
          ! units:
          call MDF_Put_Att( mf%hid, mf%varid_ap_bounds, 'units', 'Pa', status )
          IF_NOTOK_RETURN(status=1)

          ! variable:
          call MDF_Def_Var( mf%hid, 'b_bounds', MDF_DOUBLE, (/mf%dimid_nv,mf%dimid_lev/), mf%varid_b_bounds, status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_b_bounds, 'long_name', 'atmosphere hybrid sigma pressure coefficient&
               & b at layer interfaces', status )
          IF_NOTOK_RETURN(status=1)
          ! units:
          call MDF_Put_Att( mf%hid, mf%varid_b_bounds, 'units', '1', status )
          IF_NOTOK_RETURN(status=1)
        
        ! layer interfaces
        case ( 'w' )

          ! hybride layers:
          call MDF_Def_Dim( mf%hid, 'lev', levi%nlev+1, mf%dimid_lev, status )
          IF_NOTOK_RETURN(status=1)

          ! variable:
          call MDF_Def_Var( mf%hid, 'lev', MDF_DOUBLE, (/mf%dimid_lev/), mf%varid_lev, status )
          IF_NOTOK_RETURN(status=1)
          ! convert to CF name and unit:
          cf_standard_name = 'atmosphere_hybrid_sigma_pressure_coordinate'
          call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
          IF_NOTOK_RETURN(status=1)
          ! add standard_name and units attributes:
          call CF_Put_Standard_Atts( mf%hid, mf%varid_lev, trim(cf_standard_name), trim(cf_units), status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'long_name', 'pressure at layer interfaces', status )
          IF_NOTOK_RETURN(status=1)
          ! direction of increasing pressure:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'positive', 'down', status )
          IF_NOTOK_RETURN(status=1)
          ! how to compute pressure:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'forumula', 'p(n,k,j,i) = ap(k) + b(k)*ps(n,j,i)', status )
          IF_NOTOK_RETURN(status=1)
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'forumula_terms', 'ap: ap b: b ps: ps', status )
          IF_NOTOK_RETURN(status=1)
          ! generic axis name:
          call MDF_Put_Att( mf%hid, mf%varid_lev, 'axis', 'Z', status )
          IF_NOTOK_RETURN(status=1)

          ! variable:
          call MDF_Def_Var( mf%hid, 'ap', MDF_DOUBLE, (/mf%dimid_lev/), mf%varid_ap, status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_ap, 'long_name', 'atmosphere hybrid sigma pressure coefficient ap&
               & at layer interfaces', status )
          IF_NOTOK_RETURN(status=1)
          ! units:
          call MDF_Put_Att( mf%hid, mf%varid_ap, 'units', 'Pa', status )
          IF_NOTOK_RETURN(status=1)

          ! variable:
          call MDF_Def_Var( mf%hid, 'b', MDF_DOUBLE, (/mf%dimid_lev/), mf%varid_b, status )
          IF_NOTOK_RETURN(status=1)
          ! description:
          call MDF_Put_Att( mf%hid, mf%varid_b, 'long_name', 'atmosphere hybrid sigma pressure coefficient&
               & b at layer interfaces', status )
          IF_NOTOK_RETURN(status=1)
          ! units:
          call MDF_Put_Att( mf%hid, mf%varid_b, 'units', '1', status )
          IF_NOTOK_RETURN(status=1)
        
        ! other ...
        case default
        
          write (gol,'("unsupported nw : ",a)') trim(nw); call goErr
          TRACEBACK; status=1; return
          
      end select
        
      ! surface pressure:
      call MDF_Def_Var( mf%hid, 'ps', MDF_FLOAT, (/mf%dimid_lon ,mf%dimid_lat ,mf%dimid_time/), mf%varid_ps, status )
      IF_NOTOK_RETURN(status=1)
      ! convert to CF name and unit:
      cf_standard_name = 'surface_air_pressure'
      call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
      IF_NOTOK_RETURN(status=1)
      ! add standard_name and units attributes:
      call CF_Put_Standard_Atts( mf%hid, mf%varid_ps, trim(cf_standard_name), trim(cf_units), status )
      IF_NOTOK_RETURN(status=1)
      ! for post processing ...
      l = len_trim(cell_measure)
      if ( l > 0 ) then
        call MDF_Put_Att( mf%hid, mf%varid_ps, 'cell_measure', cell_measure(2:l), status )
        IF_NOTOK_RETURN(status=1)
      end if
      ! specify how computed:
      l = len_trim(cell_methods)
      if ( l > 0 ) then
        call MDF_Put_Att( mf%hid, mf%varid_ps, 'cell_methods', cell_methods(2:l), status )
        IF_NOTOK_RETURN(status=1)
      end if
      ! extra coordinates that apply to this variable:
      l = len_trim(coordinates)
      if ( l > 0 ) then
        call MDF_Put_Att( mf%hid, mf%varid_ps, 'coordinates', coordinates(2:l), status )
        IF_NOTOK_RETURN(status=1)
      end if
      
      ! extra:
      if ( nuv /= 'n' ) then
        ! surface pressure:
        call MDF_Def_Var( mf%hid, 'ps_u', MDF_FLOAT, (/mf%dimid_lonb,mf%dimid_lat ,mf%dimid_time/), mf%varid_ps_u, status )
        IF_NOTOK_RETURN(status=1)
        ! convert to CF name and unit:
        cf_standard_name = 'surface_air_pressure'
        call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
        IF_NOTOK_RETURN(status=1)
        ! add standard_name and units attributes:
        call CF_Put_Standard_Atts( mf%hid, mf%varid_ps_u, trim(cf_standard_name), trim(cf_units), status )
        IF_NOTOK_RETURN(status=1)
        
        ! surface pressure:
        call MDF_Def_Var( mf%hid, 'ps_v', MDF_FLOAT, (/mf%dimid_lon ,mf%dimid_latb,mf%dimid_time/), mf%varid_ps_v, status )
        IF_NOTOK_RETURN(status=1)
        ! convert to CF name and unit:
        cf_standard_name = 'surface_air_pressure'
        call TMM_CF_Standard_Units( trim(cf_standard_name), cf_units, status )
        IF_NOTOK_RETURN(status=1)
        ! add standard_name and units attributes:
        call CF_Put_Standard_Atts( mf%hid, mf%varid_ps_v, trim(cf_standard_name), trim(cf_units), status )
        IF_NOTOK_RETURN(status=1)
      end if
      
    end if  ! levels
    
    ! * meteo variables
    
    ! copy the paramkeys ('-aa-bb-cc-') except for the first '-' :
    l = len_trim(mf%paramkeys)
    line = mf%paramkeys(2:l)
    ! loop over all parameters:
    do iparam = 1, mf%nparam

      ! split at '-', read first part:
      call goReadFromLine( line, paramkey, status, sep='-' )
      IF_NOTOK_RETURN(status=1)

      ! define variable:
      if ( present(levi) ) then
        ! 3D field:
        if ( trim(paramkey) == 'mfu' ) then
          call MDF_Def_Var( mf%hid, trim(paramkey), MDF_FLOAT, (/mf%dimid_lonb,mf%dimid_lat ,mf%dimid_lev_u,mf%dimid_time/), &
               varid, status )
          IF_NOTOK_RETURN(status=1)
        else if ( trim(paramkey) == 'mfv' ) then
          call MDF_Def_Var( mf%hid, trim(paramkey), MDF_FLOAT, (/mf%dimid_lon ,mf%dimid_latb,mf%dimid_lev_v,mf%dimid_time/), &
               varid, status )
          IF_NOTOK_RETURN(status=1)
        else
          call MDF_Def_Var( mf%hid, trim(paramkey), MDF_FLOAT, (/mf%dimid_lon ,mf%dimid_lat ,mf%dimid_lev,mf%dimid_time/), &
               varid, status )
          IF_NOTOK_RETURN(status=1)
        end if
      else
        ! 2D field:
        if ( trim(mf%tres) == 'constant' ) then
          call MDF_Def_Var( mf%hid, trim(paramkey), MDF_FLOAT, (/mf%dimid_lon,mf%dimid_lat/), varid, status )
          IF_NOTOK_RETURN(status=1)
        else
          call MDF_Def_Var( mf%hid, trim(paramkey), MDF_FLOAT, (/mf%dimid_lon,mf%dimid_lat,mf%dimid_time/), varid, status )
          IF_NOTOK_RETURN(status=1)
        end if
      end if

      ! convert from tm5 name to standard name:
      call TMM_CF_Convert_Name( trim(paramkey), mf%cfname_param(iparam), status )
      IF_NOTOK_RETURN(status=1)
      ! get standard units:
      call TMM_CF_Standard_Units( trim(mf%cfname_param(iparam)), mf%cfunit_param(iparam), status )
      IF_NOTOK_RETURN(status=1)
      ! add standard_name and units attributes:
      call CF_Put_Standard_Atts( mf%hid, varid, trim(mf%cfname_param(iparam)), trim(mf%cfunit_param(iparam)), status )
      IF_NOTOK_RETURN(status=1)
      ! for post processing ...
      l = len_trim(cell_measure)
      if ( l > 0 ) then
        call MDF_Put_Att( mf%hid, varid, 'cell_measure', cell_measure(2:l), status )
        IF_NOTOK_RETURN(status=1)
      end if
      ! specify how computed:
      l = len_trim(cell_methods)
      if ( l > 0 ) then
        call MDF_Put_Att( mf%hid, varid, 'cell_methods', cell_methods(2:l), status )
        IF_NOTOK_RETURN(status=1)
      end if
      ! extra coordinates that apply to this variable:
      l = len_trim(coordinates)
      if ( l > 0 ) then
        call MDF_Put_Att( mf%hid, varid, 'coordinates', coordinates(2:l), status )
        IF_NOTOK_RETURN(status=1)
      end if

      ! store variable id:
      mf%varid_param(iparam) = varid
      ! store variable name:
      mf%varname_param(iparam) = trim(paramkey)
      ! no records written yet:
      mf%itrec_param(iparam) = 0

    end do   

    !
    ! end of defintion phase:
    !
    
    call MDF_EndDef( mf%hid, status )
    IF_NOTOK_RETURN(status=1)
    
    !
    ! fill time independent variables
    !
    
    ! axis:
    call MDF_Put_Var( mf%hid, mf%varid_lon, lli%lon_deg, status )
    IF_NOTOK_RETURN(status=1)
    ! axis bounds:
    call MDF_Put_Var( mf%hid, mf%varid_lon_bounds, lli%lon_bounds_deg, status )
    IF_NOTOK_RETURN(status=1)
    ! extra:
    if ( nuv /= 'n' ) then
      call MDF_Put_Var( mf%hid, mf%varid_lonb, lli%blon_deg, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! axis:
    call MDF_Put_Var( mf%hid, mf%varid_lat, lli%lat_deg, status )
    IF_NOTOK_RETURN(status=1)
    ! axis bounds:
    call MDF_Put_Var( mf%hid, mf%varid_lat_bounds, lli%lat_bounds_deg, status )
    IF_NOTOK_RETURN(status=1)
    ! extra:
    if ( nuv /= 'n' ) then
      call MDF_Put_Var( mf%hid, mf%varid_latb, lli%blat_deg, status )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! cell area ?
    if ( nuv == 'n' ) then
      ! storage:
      allocate( pat(lli%nlon,lli%nlat) )
      ! fill:
      call AreaOper( lli, pat, '=', 'm2', status )
      IF_NOTOK_RETURN(status=1)
      ! write:
      call MDF_Put_Var( mf%hid, mf%varid_cell_area, pat, status )
      IF_NOTOK_RETURN(status=1)
      ! clear:
      deallocate( pat )
    end if

    ! levels ?
    if ( present(levi) ) then
      !! tm5 variables:
      !call MDF_Put_Var( mf%hid, mf%varid_tm5_at, levi%a, status )
      !IF_NOTOK_RETURN(status=1)
      !call MDF_Put_Var( mf%hid, mf%varid_tm5_bt, levi%b, status )
      !IF_NOTOK_RETURN(status=1)
      ! where defined ?
      select case ( nw )
        !~ layer mid points:
        case ( 'n' )!, '*' )
          ! number of layers:
          !if ( present(nlev) ) then
          !  n = nlev
          !else
            n = levi%nlev
          !end if
          ! standard pressures (full levels):
          call MDF_Put_Var( mf%hid, mf%varid_lev, levi%fp0(1:n), status )
          IF_NOTOK_RETURN(status=1)
          ! full level coefficients:
          call MDF_Put_Var( mf%hid, mf%varid_ap, levi%fa(1:n), status )
          IF_NOTOK_RETURN(status=1)
          call MDF_Put_Var( mf%hid, mf%varid_b , levi%fb(1:n), status )
          IF_NOTOK_RETURN(status=1)
          ! half level coefficients:
          call MDF_Put_Var( mf%hid, mf%varid_ap_bounds, levi%fa_bounds(:,1:n), status )
          IF_NOTOK_RETURN(status=1)
          call MDF_Put_Var( mf%hid, mf%varid_b_bounds , levi%fb_bounds(:,1:n), status )
          IF_NOTOK_RETURN(status=1)
          ! extra ?
          if ( nuv /= 'n' ) then
            ! standard pressures (full levels):
            call MDF_Put_Var( mf%hid, mf%varid_lev_u, levi%fp0(1:n), status )
            IF_NOTOK_RETURN(status=1)
            ! standard pressures (full levels):
            call MDF_Put_Var( mf%hid, mf%varid_lev_v, levi%fp0(1:n), status )
            IF_NOTOK_RETURN(status=1)
          end if
        !~ layer interfaces:
        case ( 'w' )
          ! standard pressures (half levels):
          call MDF_Put_Var( mf%hid, mf%varid_lev, levi%p0, status )
          IF_NOTOK_RETURN(status=1)
          ! half level coefficients:
          call MDF_Put_Var( mf%hid, mf%varid_ap, levi%a, status )
          IF_NOTOK_RETURN(status=1)
          call MDF_Put_Var( mf%hid, mf%varid_b , levi%b, status )
          IF_NOTOK_RETURN(status=1)
        !~ other ...
        case default
          write (gol,'("unsupported nw : ",a)') trim(nw); call goErr
          TRACEBACK; status=1; return
      end select
    end if
    
    !
    ! done
    !    
    
    ! ok
    status = 0

  end subroutine Define_File
  
  
  ! ***
  
  
  subroutine WriteTimes( mf, itrec, tref, t1, t2, status )

    use GO      , only : TDate
    use GO      , only : operator(-), operator(+), operator(/), rTotal, Get, NewDate, IsAnyDate
    use MDF     , only : MDF_Put_Var

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), intent(inout)  ::  mf
    integer, intent(in)                  ::  itrec
    type(TDate), intent(in)              ::  tref, t1, t2
    integer, intent(out)                 ::  status

    ! --- const --------------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/WriteTimes'
    
    ! --- local ------------------------------
    
    type(TDate)   ::  tsince
    type(TDate)   ::  tmid
    real(8)       ::  tsec
    integer       ::  time6(6)

    ! --- begin ---------------------------------
    
    ! not for constant fields ...
    if ( trim(mf%tres) /= 'constant' ) then
  
      ! base time:
      tsince = NewDate( time6=since_time6 )

      ! write reference times:
      if ( IsAnyDate(tref) ) then
        tsec = 0.0
        time6 = 0
      else
        tsec = rTotal( tref - tsince, 'sec' )
        call Get( tref, time6=time6 )
      end if
      call MDF_Put_Var( mf%hid, mf%varid_reftime, (/tsec/), status, start=(/itrec/), count=(/1/) )
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Var( mf%hid, mf%varid_reftimevalues, time6, status, start=(/1,itrec/), count=(/6,1/) )
      IF_NOTOK_RETURN(status=1)

      ! write mid time:
      if ( IsAnyDate(t1) .or. IsAnyDate(t2) ) then
        tsec = 0.0
        time6 = 0
      else
        tmid = t1 + (t2-t1)/2
        tsec = rTotal( tmid - tsince, 'sec' )
        call Get( tmid, time6=time6 )
      end if
      call MDF_Put_Var( mf%hid, mf%varid_time, (/tsec/), status, start=(/itrec/), count=(/1/) )
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Var( mf%hid, mf%varid_timevalues, time6, status, start=(/1,itrec/), count=(/6,1/) )
      IF_NOTOK_RETURN(status=1)

      ! start time:
      if ( IsAnyDate(t1) ) then
        tsec = 0.0
        time6 = 0
      else
        tsec = rTotal( t1 - tsince, 'sec' )
        call Get( t1, time6=time6 )
      end if
      call MDF_Put_Var( mf%hid, mf%varid_time_bounds, (/tsec/), status, start=(/1,itrec/), count=(/1,1/) )
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Var( mf%hid, mf%varid_timevalues_bounds, time6, status, start=(/1,1,itrec/), count=(/1,6,1/) )
      IF_NOTOK_RETURN(status=1)

      ! end time:
      if ( IsAnyDate(t2) ) then
        tsec = 0.0
        time6 = 0
      else
        tsec = rTotal( t2 - tsince, 'sec' )
        call Get( t2, time6=time6 )
      end if
      call MDF_Put_Var( mf%hid, mf%varid_time_bounds, (/tsec/), status, start=(/2,itrec/), count=(/1,1/) )
      IF_NOTOK_RETURN(status=1)
      call MDF_Put_Var( mf%hid, mf%varid_timevalues_bounds, time6, status, start=(/2,1,itrec/), count=(/1,6,1/) )
      IF_NOTOK_RETURN(status=1)
      
    end if   ! not const
    
    ! ok
    status = 0
    
  end subroutine WriteTimes
  
    
  ! ***
  
  
  subroutine mf_WriteRecord_2d( mf, tmi, paramkey, unit, tref, t1, t2, &
                                lli, nuv, ll, status )

    use GO      , only : TDate
    use Grid    , only : TllGridInfo
    use MDF     , only : MDF_Create, MDF_Close, MDF_REPLACE
    use MDF     , only : MDF_Put_Var
    use tmm_info, only : TMeteoInfo
    use TMM_CF  , only : TMM_CF_Convert_Units

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), 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 ------------------------------
    
    integer       ::  iparam
    integer       ::  i
    integer       ::  itrec
    type(TDate)   ::  tsince
    type(TDate)   ::  tmid
    real(8)       ::  tsec
    integer       ::  time6(6)
    real          ::  ufac

    ! --- 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 MDF_Create( trim(mf%fname), output_type, MDF_REPLACE, mf%hid, status )
      IF_NOTOK_RETURN(status=1)
      ! write file header:
      call Define_File( mf, lli, nuv, 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.
    endif

    ! search parameter:
    iparam = -1
    do i = 1, mf%nparam
      if ( trim(paramkey) == trim(mf%varname_param(i)) ) then
        iparam = i
        exit
      end if
    end do
    if ( iparam < 0 ) then
      write (gol,'("could not find parameter `",a,"` in list : ",a)') trim(paramkey), trim(mf%paramkeys); call goErr
      TRACEBACK; status=1; return
    end if
    
    ! next record:
    mf%itrec_param(iparam) = mf%itrec_param(iparam) + 1
    ! short:
    itrec = mf%itrec_param(iparam)

    ! update time fields:
    call WriteTimes( mf, itrec, tref, t1, t2, status )
    IF_NOTOK_RETURN(status=1)
    
    ! conversion factor:
    call TMM_CF_Convert_Units( trim(unit), trim(mf%cfunit_param(iparam)), ufac, status )
    IF_NOTOK_RETURN(status=1)
    !! info ...
    !if ( ufac /= 1.0 ) then
    !  write (gol,'("      convert `",a,"` from `",a,"` to `",a,"` with factor ",f8.2," ; write range ",2f8.2)') &
    !           trim(paramkey), trim(unit), trim(mf%cfunit_param(iparam)), ufac, minval(ll*ufac), maxval(ll*ufac); call goPr
    !end if
 
    ! write data:
    if ( trim(mf%tres) == 'constant' ) then
      call MDF_Put_Var( mf%hid, mf%varid_param(iparam), ll*ufac, status )
      IF_NOTOK_RETURN(status=1)
    else
      call MDF_Put_Var( mf%hid, mf%varid_param(iparam), ll*ufac, status, &
                            start=(/1,1,itrec/), count=(/size(ll,1),size(ll,2),1/) )
      IF_NOTOK_RETURN(status=1)
    end if
 
    ! 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%nparam ) then
      ! close file:
      call MDF_Close( mf%hid, status )
      IF_NOTOK_RETURN(status=1)
      ! rewrite status file:
      call WriteStatus( mf, 'completed', status )
      IF_NOTOK_RETURN(status=1)
      ! reset flag:
      mf%output_initialised = .false.
    end if

    ! 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 MDF     , only : MDF_Create, MDF_Close, MDF_REPLACE
    use MDF     , only : MDF_Put_Var
    use tmm_info, only : TMeteoInfo
    use TMM_CF  , only : TMM_CF_Convert_Units

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), 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 ------------------------------
    
    integer             ::  iparam
    integer             ::  i
    integer             ::  itrec
    type(TDate)         ::  tsince
    type(TDate)         ::  tmid
    real(8)             ::  tsec
    integer             ::  time6(6)
    real                ::  ufac
    real, allocatable   ::  llX(:,:,:)
    character(len=1)    ::  nwX

    ! --- 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
    
    ! convection fields provided on lowest layers only ...
    if ( nw == '*' ) then
      ! extend to all layers:
      allocate( llX(size(ll,1),size(ll,2),levi%nlev) )
      ! pad with zeros:
      llX = 0.0
      llX(:,:,1:size(ll,3)) = ll
      ! new description:
      nwX = 'n'
    else
      ! just copy:
      allocate( llX(size(ll,1),size(ll,2),size(ll,3)) )
      llX = ll
      nwX = nw
    end if

    ! new or existing ?
    if ( .not. mf%output_initialised ) then
      ! open new file, destroy old:
      call MDF_Create( trim(mf%fname), output_type, MDF_REPLACE, mf%hid, status )
      IF_NOTOK_RETURN(status=1)
      ! write file header:
      call Define_File( mf, lli, nuv, status, levi=levi, nw=nwX )!, nlev=nlev )
      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.
    endif
      
    ! search parameter:
    iparam = -1
    do i = 1, mf%nparam
      if ( trim(paramkey) == trim(mf%varname_param(i)) ) then
        iparam = i
        exit
      end if
    end do
    if ( iparam < 0 ) then
      write (gol,'("could not find parameter `",a,"` in list : ",a)') trim(paramkey), trim(mf%paramkeys); call goErr
      TRACEBACK; status=1; return
    end if
    
    ! next record:
    mf%itrec_param(iparam) = mf%itrec_param(iparam) + 1
    ! short:
    itrec = mf%itrec_param(iparam)

    ! update time fields:
    call WriteTimes( mf, itrec, tref, t1, t2, status )
    IF_NOTOK_RETURN(status=1)
    
    ! conversion factor:
    call TMM_CF_Convert_Units( trim(unit), trim(mf%cfunit_param(iparam)), ufac, status )
    IF_NOTOK_RETURN(status=1)
    !! info ...
    !if ( ufac /= 1.0 ) then
    !  write (gol,'("      convert `",a,"` from `",a,"` to `",a,"` with factor ",f8.2," ; write range ",2f8.2)') &
    !           trim(paramkey), trim(unit), trim(mf%cfunit_param(iparam)), ufac, minval(llX*ufac), maxval(llX*ufac); call goPr
    !end if
    ! convert:
    llX = llX * ufac
    ! write data:
    call MDF_Put_Var( mf%hid, mf%varid_param(iparam), llX, status, &
                          start=(/1,1,1,itrec/), count=(/size(llX,1),size(llX,2),size(llX,3),1/) )
    IF_NOTOK_RETURN(status=1)
    
    ! clear:
    deallocate( llX )

    ! write surface pressure:
    if ( nuv == 'u' ) then
      call MDF_Put_Var( mf%hid, mf%varid_ps_u, ps, status, &
                            start=(/1,1,itrec/), count=(/size(ps,1),size(ps,2),1/) )
      IF_NOTOK_RETURN(status=1)
    else if ( nuv == 'v' ) then
      call MDF_Put_Var( mf%hid, mf%varid_ps_v, ps, status, &
                            start=(/1,1,itrec/), count=(/size(ps,1),size(ps,2),1/) )
      IF_NOTOK_RETURN(status=1)
    else
      call MDF_Put_Var( mf%hid, mf%varid_ps, ps, status, &
                            start=(/1,1,itrec/), count=(/size(ps,1),size(ps,2),1/) )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! 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%nparam ) then
      ! close file:
      call MDF_Close( mf%hid, status )
      IF_NOTOK_RETURN(status=1)
      ! rewrite status file:
      call WriteStatus( mf, 'completed', status )
      IF_NOTOK_RETURN(status=1)
      ! reset flag:
      mf%output_initialised = .false.
    end if

    ! ok
    status = 0

  end subroutine mf_WriteRecord_3d
  

  ! ***
  
  
  subroutine WriteStatus( mf, msg, status )

    ! --- in/out -------------------------------
    
    type(TMeteoFile_tm5_nc), 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
  

end module tmm_mf_tm5_nc