ecearth3/sources/tm5mp/proj/co2/user_output_pdump__co2.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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
#define IF_NOTOK_MDF(action) if (status/=0) then; TRACEBACK; action; call MDF_CLose(fid,status); status=1; return; end if
!
#include "tm5.inc"
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      USER_OUTPUT_PDUMP
!
! !DESCRIPTION: 
!
! Module to deal with time-series output. Output are in NetCDF-4 and use CF
! conventions. The following output are available:
! 
!     - one file with grid definition
!     - one file with time series of some met fields (pressure, temperature, winds, ...)
!     - one or more files with time series of some tracers
!     - one or two files with Local Time output for some tracers
!     - one file with time series of wet and dry depositions
!     - one file with time series of deposition velocity
!
! Activation, tracers to account for, time step of the series, are set in the
! rcfile, following this template :
!     
!
! SAMPLE RCFILE
!
!    output.pdump                        :  T
!    output.pdump.dataset.author         :  John Doe
!    output.pdump.dataset.institution    :  MyFirm, Anytown, USA
!    output.pdump.dataset.version        :  GEMS GRG; era2003 simulation
!    output.pdump.fname.model            :  TM5
!    output.pdump.fname.expid            :  V2
!    output.pdump.fname.grid.300x200     :  3x2   ! short name, required if there is zoom regions
!    output.pdump.fname.grid.100x100     :  1x1
!
!    output.pdump.griddef.apply          :  T
!
!    output.pdump.tp.apply               :  T
!    output.pdump.tp.dhour               :  1
!
!    output.pdump.vmr.n                  :  3
!
!    output.pdump.vmr.001.apply          :  T
!    output.pdump.vmr.001.fname          :  vmr1
!    output.pdump.vmr.001.dhour          :  3
!    output.pdump.vmr.001.tracers        :  CO2
!
!    output.pdump.lt.apply               :  T
!    output.pdump.lt.tracers             :  CO2
!    output.pdump.lt.localtime           :  2
!    
!    output.pdump.lt2.apply              :  F
!    output.pdump.lt2.tracers            :  CO2
!    output.pdump.lt2.localtime          :  12
! 
!    output.pdump.depositions.apply      :  F
!    output.pdump.depositions.dhour      :  3
!    output.pdump.depositions.tracers    :  CO2
!
!    output.pdump.depvels.apply          :  F
!    output.pdump.depvels.dhour          :  3
!    output.pdump.depvels.tracers        :  CO2
!
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_PDUMP
  !
  ! !USES:
  !
  use GO,           only : gol, goPr, goErr, goLabel
  use GO,           only : TDate
  use dims,         only : nregions, idatee, idatei, okdebug
  use chem_param,   only : ntrace
  
  USE MDF
  USE TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo
  
  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  ::  Output_PDUMP_Init
  public  ::  Output_PDUMP_Step
  public  ::  Output_PDUMP_Done
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter :: mname            = 'user_output_pdump'
  character(len=*), parameter :: outfileversnr    = '0.1'
  integer,          parameter :: time_reftime6(6) = (/1950,01,01,00,00,00/)          ! reference time
  character(len=*), parameter :: time_units       = 'days since 1950-01-01 00:00:00' 
  !
  !
  type TPdumpFile_GridDef
    integer             ::  trec
    integer             ::  ncid
    integer             ::  dimid_scalar, dimid_lon, dimid_lat, dimid_lev, dimid_levi
    integer             ::  varid_lon, varid_lat, varid_time, varid_date
    integer             ::  varid_gridbox_area
    integer             ::  varid_a, varid_b
    integer             ::  varid_a_bnds, varid_b_bnds
    integer             ::  varid_p0
    !integer             ::  varid_ps
    !integer             ::  varid_geo_height
  end type TPdumpFile_GridDef

  type TPdumpFile_TP
    integer             ::  trec
    integer             ::  dhour
    integer             ::  ncid
    integer             ::  dimid_lon, dimid_lat, dimid_lev, dimid_time, dimid_datelen
    integer             ::  varid_lon, varid_lat, varid_lev, varid_time, varid_date
    integer             ::  varid_ps
    integer             ::  varid_surface_temp
    integer             ::  varid_orog
    integer             ::  varid_geop
    integer             ::  varid_pressure
    integer             ::  varid_temp
    integer             ::  varid_humid
    integer             ::  varid_u, varid_v, varid_w
    real, allocatable   ::  data3d(:,:,:,:,:)
    real, allocatable   ::  data2d(:,:,:,:)
    real, allocatable   ::  time(:)
    real, allocatable   ::  date(:,:)
  end type TPdumpFile_TP

  type TPdumpFile_VMR
     integer             ::  trec, n_rec
     logical             :: apply
    integer             ::  dhour
    character(len=256)  ::  tracer_names
    integer             ::  ncid
    integer             ::  dimid_lon, dimid_lat, dimid_lev, dimid_levi, dimid_time, dimid_datelen
    integer             ::  varid_lon, varid_lat, varid_lev, varid_time, varid_date
    integer             ::  varid_temp    
    integer             ::  varid_ps
    integer             ::  varid_a_bnds, varid_b_bnds
    integer             ::  ntr
    integer             ::  itr(ntrace)
    character(len=8)    ::  name_tr(ntrace)
    integer             ::  varid_tr(ntrace)
    character(len=4)    ::  varid_type(ntrace)
    real, allocatable   ::  data3d(:,:,:,:,:)
    real, allocatable   ::  data3d_t(:,:,:,:)
    real, allocatable   ::  sp(:,:,:)
    real, allocatable   ::  time(:)
    real, allocatable   ::  date(:,:)
  end type TPdumpFile_VMR

  type TPdumpFile_LT
    integer             ::  trec
    character(len=256)  ::  tracer_names
    integer             ::  ncid
    integer             ::  local_time
    integer             ::  dimid_lon, dimid_lat, dimid_lev, dimid_time, dimid_datelen
    integer             ::  varid_lon, varid_lat, varid_lev, varid_time, varid_date
    integer             ::  varid_ps
    integer             ::  ntr
    integer             ::  itr(ntrace)
    character(len=8)    ::  name_tr(ntrace)
    integer             ::  varid_tr(ntrace)
    real,allocatable    ::  accu(:,:,:,:)
    real,allocatable    ::  naccu(:,:)
    real,allocatable    ::  p_accu(:,:)
    real,allocatable    ::  np_accu(:)
  end type TPdumpFile_LT

  type TPdumpFile_DEPS
    integer             ::  trec
    integer             ::  dhour
    character(len=256)  ::  tracer_names
    integer             ::  ncid
    integer             ::  dimid_lon, dimid_lat, dimid_time, dimid_datelen
    integer             ::  varid_lon, varid_lat, varid_time, varid_date, varid_accum
    integer             ::  ntr
    integer             ::  itr(ntrace)
    character(len=8)    ::  name_tr(ntrace)
    integer             ::  varid_ddep(ntrace)
    real, pointer       ::  ddep_budget(:,:,:)
    logical             ::  with_wdep(ntrace)
    integer             ::  varid_wdep(ntrace)
    real, pointer       ::  wdep_budget(:,:,:)
    type(TDate)         ::  t0_budget
    real, allocatable   ::  data2d_dry(:,:,:,:)
    real, allocatable   ::  data2d_wet(:,:,:,:)
    real, allocatable   ::  time(:), dt(:)
    real, allocatable   ::  date(:,:)
  end type TPdumpFile_DEPS

  type TPdumpFile_DEPV
    integer             ::  trec
    integer             ::  dhour
    character(len=256)  ::  tracer_names
    integer             ::  ncid
    integer             ::  dimid_lon, dimid_lat, dimid_time, dimid_datelen
    integer             ::  varid_lon, varid_lat, varid_time, varid_date
    integer             ::  ntr
    integer             ::  itr(ntrace)
    character(len=8)    ::  name_tr(ntrace)
    integer             ::  varid_tr(ntrace)
    real, allocatable   ::  data2d(:,:,:,:)
    real, allocatable   ::  time(:)
    real, allocatable   ::  date(:,:)
  end type TPdumpFile_DEPV


  ! --- var -----------------------------

  integer                         ::  fid ! file id for IF_NOTOK_MDF macro
  integer                         ::  access_mode ! netcdf-4 access mode
  integer                         ::  curr_day(nregions,3)

  logical                         ::  firstday
  logical                         ::  lastday  ! it is last day and not a full day (ie day does not end at 00 of next day) 

  character(len=32)               ::  fname_model
  character(len=8)                ::  fname_expid, meteo_class
  character(len=32)               ::  fname_grid(nregions)
  character(len=256)              ::  dataset_author, institution, dataset_version

  logical, save                                ::  griddef_apply
  type(TPdumpFile_GridDef), save               ::  RF_GridDef(nregions)

  logical, save                                ::  tp_apply
  integer                                      ::  tp_dhour, n_tp_rec
  type(TPdumpFile_TP), save                    ::  RF_TP(nregions)

  integer, save                                ::  nvmr
  logical, allocatable                         ::  vmr_apply(:)
  real, allocatable                            ::  vmr_sregbord(:,:)
  character(len=16), allocatable               ::  vmr_fname(:)
  integer, allocatable                         ::  vmr_dhour(:)
  character(len=256), allocatable              ::  vmr_tracer_names(:)
  type(TPdumpFile_VMR), allocatable, save      ::  RF_VMR(:,:)

  logical, save                                ::  lt_apply
  character(len=16)                            ::  lt_fname
  character(len=256)                           ::  lt_tracer_names
  integer                                      ::  lt_localtime
  type(TPdumpFile_LT), save                    ::  RF_LT(nregions)

  logical, save                                ::  lt2_apply
  character(len=16)                            ::  lt2_fname
  character(len=256)                           ::  lt2_tracer_names
  integer                                      ::  lt2_localtime
  type(TPdumpFile_LT), save                    ::  RF_LT2(nregions)

  logical, save                                ::  deps_apply
  character(len=16)                            ::  deps_fname
  integer                                      ::  deps_dhour, n_deps_rec
  character(len=256)                           ::  deps_tracer_names
  type(TPdumpFile_DEPS), save                  ::  RF_DEPS(nregions)

  logical, save                                ::  depv_apply
  character(len=16)                            ::  depv_fname
  integer                                      ::  depv_dhour, n_depv_rec
  character(len=256)                           ::  depv_tracer_names
  type(TPdumpFile_DEPV), save                  ::  RF_DEPV(nregions)
  !
  ! !REVISION HISTORY:
  !
  !    1 Oct 2010 - Achim Strunk - revised older RETRO ouptut :
  !                                  add 2nd local time, regional output
  !   10 Jul 2012 - Ph. Le Sager - switch from pnetcdf to netcdf4_par (through
  !                                  MDF); get rid of the with_tendencies code.
  !   12 Nov 2012 - Ph. Le Sager - adapted for lon-lat MPI decomposition.
  !                              - get rid of unlimited dimensions so we can
  !                                  write in collective mode.
  !                              - store series to write them only at end-of-day
  !                                  to speed-up code
  !   10 Oct 2013 - Ph. Le Sager - fixed GET_N_TIME_RECORDS and several 'init'
  !                                 and write' routines.
  !   14 Apr 2014 - Ph. Le Sager + JEW - tropomi add-ons in VMR: Temperature,
  !                                      As, Bs, better CF
  !
  ! !REMARKS:
  ! 
  !   (1) Initially called RETRO output for GEMS GRG, the module has been adapted
  !        for CLIMAQS project and renamed PDUMP.
  !   (2) Previous remarks "as is":
  !        - longitudes from [0,360] ?
  !          this is imposible for zoom area's such as for the heatwave
  !        - levels from surface to top
  !        - time from 1950-01-01 00:00
  !   (3) This is supposed to work with netcdf4_parallel. You cannot use
  !         MPI with a non-parallel version of netcdf4 here.
  !   (4) The parallel writing is done in COLLECTIVE mode, but remain
  !         expensive on some system. Possible optimization : output one file
  !         per month (chunk/leg), and/or per field, and/or per processor. 
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    GET_N_TIME_RECORDS
  !
  ! !DESCRIPTION: return number of time steps for a daily timeseries file
  !\\
  !\\
  ! !INTERFACE:
  !
  FUNCTION GET_N_TIME_RECORDS( date, dh, isDEPS, mess )
    !
    ! !USES:
    !
    USE GO , only : TDate, NewDate, rTotal, operator(-)
    !
    ! !RETURN VALUE:
    !
    integer :: get_n_time_records
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: date(6)  ! 1st time step of the day (date[4] should be 1 unless 1st day)
    integer, intent(in)  :: dh       ! time step for timeseries (should divide 24)
    logical,          optional, intent(in) :: isDEPS ! to differentiate b/w DEPS and others
    character(len=*), optional, intent(in) :: mess   ! message (if okdebug)
    !
    ! !REVISION HISTORY: 
    !    9 Nov 2012 - Ph. Le Sager - v0
    !    9 Oct 2013 - Ph. Le Sager - fix to work with default "output.after.step: v"
    !
    ! !REMARKS:
    !  - assumed dynamic timestep of 1h. It may be smaller at some steps (CFL
    !     violation) but that's ok, since it is never zero and we are looking
    !     at "integer-hours" for output. Only issue is if dynamic timestep is
    !     LARGER than timestep of timeseries: deemed too exotic, but we could
    !     add a test in init to avoid that.
    !
    ! !TODO:
    !  - need to check if Oct 2013 fix works with runs that do NOT stop at
    !     midnight (labelled "lastday" cases in the code). Deemed exotic for now.
    !  - check if anything changes with other possible values of "output.after.step"
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer     :: is, ie, delta, dynstep
    logical     :: deps
    type(TDate) :: t, t0
    real        :: time

    ! Type of record (standard=vmr, tp, depv) or special (deps)
    deps=.false.
    if (present(isDEPS)) deps=isDEPS

    ! Start index
    delta=date(4)
    if (deps) delta=date(4)+1  ! one DYNAMIC time step done to output something

    if (modulo(delta,dh)==0) then
       is=delta/dh
    else
       is=(delta+dh)/dh
    end if
    
    ! End index
    ie = 23 / dh    ! 23 = 24 - dynamic time step

    if (deps) then  ! there will be an extra step if run goes further than midnight
       
       t0   = NewDate( time6=date   )
       t    = NewDate( time6=idatee )
       time = rTotal( t - t0, 'day' )

       if (time > 1) ie=24/dh
    end if

    ! Case of last day stopping before midnite. (Need testing again - see
    ! !TODO. Probaly alright for VMR/TP/DEPV, but check required for DEPS)
    if (lastday) ie=idatee(4)/dh

    ! length
    get_n_time_records = ie-is+1

    if(okdebug)then
       if (present(mess))then
          write(gol,*) 'GET_N_TIME_RECORDS -'//trim(mess); call goPr
       end if
       
       write(gol,*) "GET_N_TIME_RECORDS - is, ie, deps, firstday, lastday, get_n_time_records:"; call goPr
       write(gol,*) "GET_N_TIME_RECORDS - ", is, ie, deps, firstday, lastday, get_n_time_records; call goPr
    end if
    
    return 

  END FUNCTION GET_N_TIME_RECORDS
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_PDUMP_INIT
  !
  ! !DESCRIPTION: reads rc file keys relevant for pdump
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OUTPUT_PDUMP_INIT( rcF, dhour_min, status )
    !
    ! !USES:
    !
    use GO,        only : TrcFile, ReadRc
    use MeteoData, only : lli, set
    use MeteoData, only : sp_dat,  oro_dat, temper_dat, humid_dat, pu_dat, pv_dat
    use MeteoData, only : mfw_dat, gph_dat, t2m_dat
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TrcFile), intent(inout)     ::  rcF
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  dhour_min
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - upgrade from RETRO to PDUMP
    !    8 Nov 2012 - Ph. Le Sager - added access mode switch
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer                   ::  region
    character(len=64)         ::  key
    character(len=3)          ::  nr
    integer                   ::  ivmr

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

    call goLabel(rname)

#ifdef MPI
#ifdef with_netcdf4_par       
    access_mode = MDF_COLLECTIVE
#else
    write(gol,'("Time Series output (PDUMP) requires netcdf4 with parallel access enabled")') ; call goErr
    TRACEBACK
    status=1; return
#endif
#else
    access_mode = MDF_INDEPENDENT
#endif

    ! which day
    firstday = .true.
    lastday  = .true.
    
    if (any(idatei(1:3)/=idatee(1:3))) lastday=.false.  ! i.e. at least one full day

    ! dataset keys:
    call ReadRc( rcF, 'output.pdump.dataset.author'     , dataset_author     , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.dataset.institution', institution        , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.dataset.version'    , dataset_version    , status )
    IF_NOTOK_RETURN(status=1)

    ! filename keys:
    call ReadRc( rcF, 'output.pdump.fname.model', fname_model, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.fname.expid', fname_expid, status )
    IF_NOTOK_RETURN(status=1)

    ! prefix grid name in case of zooming regions:
    if ( nregions > 1 ) then
      ! loop over regions:
      do region = 1, nregions
        ! short grid name from rcfile:
        call ReadRc( rcF, 'output.pdump.fname.grid.'//trim(lli(region)%name), key, status )
        IF_NOTOK_RETURN(status=1)
        ! fill grid extenstion to file names:
        fname_grid(region) = '-'//trim(key)
      end do
    else
      ! empty
      fname_grid = ''
    end if

    ! griddef file ?
    call ReadRc( rcF, 'output.pdump.griddef.apply', griddef_apply, status )
    IF_NOTOK_RETURN(status=1)

    ! temperature, pressure, etc file ?
    call ReadRc( rcF, 'output.pdump.tp.apply', tp_apply, status )
    IF_NOTOK_RETURN(status=1)

    if (tp_apply) then 
       ! ensure that required meteo is loaded
       do region=1,nregions
          call Set(     sp_dat(region), status, used=.true. )
          call Set(    oro_dat(region), status, used=.true. )
          call Set( temper_dat(region), status, used=.true. )
          call Set(    t2m_dat(region), status, used=.true. )
          call Set(  humid_dat(region), status, used=.true. )
          call Set(     pu_dat(region), status, used=.true. )
          call Set(     pv_dat(region), status, used=.true. )
          call Set(    mfw_dat(region), status, used=.true. )
          call Set(    gph_dat(region), status, used=.true. )   ! used to compute vertical wind
       end do
    end if

    ! time resolution (1 hour by default)
    call ReadRc( rcF, 'output.pdump.tp.dhour', tp_dhour, status, default=1 )
    IF_ERROR_RETURN(status=1)

    ! VMR files
    call ReadRc( rcF, 'output.pdump.vmr.n', nvmr, status )    ! number of vmr files to be written
    IF_NOTOK_RETURN(status=1)

    if ( nvmr < 0 ) then
       write (gol,'("strange specification of number of vmr files : ",i6)') nvmr; call goErr
       TRACEBACK; status=1; return
    end if

    ! meteo
    call ReadRc( rcF, 'my.meteo.class', meteo_class, status, default='unknown' )
    IF_ERROR_RETURN(status=1)
    
   
    ! write any vmr files ?
    if ( nvmr > 0 ) then
       ! storage:
       allocate( vmr_apply(nvmr)        ) ; vmr_apply = .false.
       allocate( vmr_fname(nvmr)        ) ; vmr_fname = ''
       allocate( vmr_dhour(nvmr)        ) ; vmr_dhour = -1
       allocate( vmr_tracer_names(nvmr) ) ; vmr_tracer_names = '' 
       allocate( vmr_sregbord(nvmr,4)   ) ; vmr_sregbord = -999.9
       allocate( RF_VMR(nregions,nvmr)  )
       
#ifdef tropomi
       do region=1,nregions
          call Set( temper_dat(region), status, used=.true. )
       end do
#endif

       ! loop over vmr files:
       do ivmr = 1, nvmr
          ! number in rc keys:
          write (nr,'(i3.3)') ivmr
          ! apply ?
          call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.apply', vmr_apply(ivmr), status )
          IF_NOTOK_RETURN(status=1)

          rf_vmr(:,ivmr)%apply = vmr_apply(ivmr)
          
          ! proceed ?
          if ( vmr_apply(ivmr) ) then
             ! first part of filename:
             call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.fname', vmr_fname(ivmr), status )
             IF_NOTOK_RETURN(status=1)
             ! time resolution:
             call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.dhour', vmr_dhour(ivmr), status )
             IF_NOTOK_RETURN(status=1)
             ! tracers to be written:
             call ReadRc( rcF, 'output.pdump.vmr.'//nr//'.tracers', vmr_tracer_names(ivmr), status )
             IF_NOTOK_RETURN(status=1)

          end if  ! apply ?
       end do  ! vmr numbers

       ! required meteo
       if (any(vmr_apply)) then
          do region=1,nregions
             call Set(     sp_dat(region), status, used=.true. )
          end do
       end if
       
    end if  ! nvmr > 0

    ! ---------------------
    ! local time:
    ! ---------------------

    ! file 1 
    lt_fname = 'lt'
    call ReadRc( rcF, 'output.pdump.lt.apply', lt_apply, status )
    IF_NOTOK_RETURN(status=1)
    if ( lt_apply ) then
      call ReadRc( rcF, 'output.pdump.lt.tracers', lt_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.pdump.lt.localtime', lt_localtime, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! file 2
    lt2_fname = 'lt2'
    call ReadRc( rcF, 'output.pdump.lt2.apply', lt2_apply, status )
    IF_NOTOK_RETURN(status=1)
    if ( lt2_apply ) then
      call ReadRc( rcF, 'output.pdump.lt2.tracers', lt2_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.pdump.lt2.localtime', lt2_localtime, status )
      IF_NOTOK_RETURN(status=1)
    end if

    if (lt_apply .or. lt2_apply) then
       do region=1,nregions
          call Set( sp_dat(region), status, used=.true. )
       end do
    end if

    ! ---------------------
    ! deposition fluxes:
    ! ---------------------
    deps_fname = 'depositions'
    call ReadRc( rcF, 'output.pdump.depositions.apply', deps_apply, status )
    IF_NOTOK_RETURN(status=1)
    if ( deps_apply ) then
      call ReadRc( rcF, 'output.pdump.depositions.dhour', deps_dhour, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.pdump.depositions.tracers', deps_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! ---------------------
    ! deposition velocities:
    ! ---------------------
    depv_fname = 'depvels'
    call ReadRc( rcF, 'output.pdump.depvels.apply', depv_apply, status )
    IF_NOTOK_RETURN(status=1)
    if ( depv_apply ) then
      call ReadRc( rcF, 'output.pdump.depvels.dhour', depv_dhour, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.pdump.depvels.tracers', depv_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! no files open yet
    curr_day = -1

    ! lowest time frequency is 1 hour
    dhour_min = 1

    call goLabel()

    ! ok
    status = 0

  END SUBROUTINE OUTPUT_PDUMP_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_PDUMP_STEP
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OUTPUT_PDUMP_STEP( region, idate_f, status )
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !
    ! !REMARKS:
    !  (1) called every hour.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    call goLabel(rname)
    
    !----------------------
    ! close if necessary
    !----------------------
    ! if a file is open, and it is a new day
    if ( all(curr_day(region,:) > 0) .and. any(idate_f(1:3) /= curr_day(region,:)) ) then

       ! write in previous day file end-of-interval data
       call PDUMP_Files_Write2( region, idate_f, status )
       IF_NOTOK_RETURN(status=1)

       ! close all
       call PDUMP_Files_Close( region, status )
       IF_NOTOK_RETURN(status=1)

       ! no files open ...
       curr_day(region,:) = -1

       firstday = .false.

    end if

    !----------------------
    ! open if necessary
    !----------------------
    if ( any(curr_day(region,:) < 0) ) then

       if (all(idate_f(1:3)==idatee(1:3))) lastday=.true.  ! means last day is not a full day

       write(gol,*) "U_O_Pdump open [idate_f, last day] = ", idate_f, lastday; call goPr

       call PDUMP_Files_Open( region, idate_f, status )
       IF_NOTOK_RETURN(status=1)

       ! store date of current day
       curr_day(region,:) = idate_f(1:3)

    end if

    !----------------------
    ! write
    !----------------------
    call PDUMP_Files_Write( region, idate_f, status )
    IF_NOTOK_RETURN(status=1)

    ! if not midnight, write end-of-interval data
    if ( any(idate_f(4:6) > 0) ) then
       call PDUMP_Files_Write2( region, idate_f, status )
       IF_NOTOK_RETURN(status=1)
    end if

    !----------------------
    ! done
    !----------------------
    call goLabel()

    status = 0

  END SUBROUTINE OUTPUT_PDUMP_STEP
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_PDUMP_DONE
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE OUTPUT_PDUMP_DONE( status )
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !   31 Aug 2012 - P. Le Sager - reverse order in which regions are dealt with (MDF issue)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Output_PDUMP_Done'
    integer                     ::  region

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

    ! close files:
    do region = nregions, 1, -1

      call PDUMP_Files_Close( region, status )
      IF_NOTOK_RETURN(status=1)

    end do

    ! clear:
    if ( nvmr > 0 ) then
      deallocate( vmr_apply        )
      deallocate( vmr_fname        )
      deallocate( vmr_dhour        )
      deallocate( vmr_tracer_names )
      deallocate( vmr_sregbord     )
      deallocate( RF_VMR           )
    end if

    ! ok
    status = 0

  END SUBROUTINE OUTPUT_PDUMP_DONE
  !EOC



  ! ********************************************************************
  ! ***
  ! *** open/write/close pdump files
  ! ***
  ! ********************************************************************


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PDUMP_FILES_OPEN
  !
  ! !DESCRIPTION: call init method of each output file.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine PDUMP_Files_Open( region, idate_f, status )
    !
    ! !USES:
    !
    use global_data, only : outdir
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer            ::  ivmr

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

    ! grid definition:
    if ( griddef_apply ) then
      call RF_GridDef_Init( RF_GridDef(region), outdir, fname_model, fname_expid, region, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! dynamics:
    if ( tp_apply ) then
      call RF_TP_Init     ( RF_TP(region)     , outdir, fname_model, fname_expid, &
                            region, idate_f, tp_dhour, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! tracer concentrations:
    do ivmr = 1, nvmr
      if ( .not. vmr_apply(ivmr) ) cycle
      call RF_VMR_Init( RF_VMR(region,ivmr),  outdir, fname_model, fname_expid, &
                        vmr_fname(ivmr),  region, idate_f,                      &
                        vmr_dhour(ivmr),  vmr_tracer_names(ivmr),   status )
      IF_NOTOK_RETURN(status=1)

      vmr_apply(ivmr) = rf_vmr(region,ivmr)%apply
    end do


    ! lt output:
    if ( lt_apply ) then
      call RF_LT_Init( RF_LT(region), outdir, fname_model, fname_expid, &
                       lt_fname, region, idate_f,                       &
                       lt_localtime, lt_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
    end if
    if ( lt2_apply ) then
      call RF_LT_Init( RF_LT2(region), outdir, fname_model, fname_expid, &
                       lt2_fname, region, idate_f,                       &
                       lt2_localtime, lt2_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! deposition fluxes:
!     if ( deps_apply ) then
!       call RF_DEPS_Init( RF_DEPS(region), outdir, fname_model, fname_expid, &
!                          deps_fname, region, idate_f, &
!                          deps_dhour, deps_tracer_names, status )
!       IF_NOTOK_RETURN(status=1)
!     end if

!     ! deposition velocities:
!     if ( depv_apply ) then
!       call RF_DEPV_Init( RF_DEPV(region), outdir, fname_model, fname_expid, &
!                          depv_fname, region, idate_f, &
!                          depv_dhour, depv_tracer_names, status )
!       IF_NOTOK_RETURN(status=1)
!     end if


    ! ok
    status = 0

  END SUBROUTINE PDUMP_FILES_OPEN
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PDUMP_FILES_WRITE
  !
  ! !DESCRIPTION: call write method for each output file.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PDUMP_FILES_WRITE( region, idate_f, status )
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer            ::  ivmr

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

    ! grid definition:
    if ( griddef_apply ) then
      call RF_GridDef_Write( RF_GridDef(region), region, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! dynamics:
    if ( tp_apply ) then
      call RF_TP_Write( RF_TP(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! tracer fields:
    do ivmr = 1, nvmr
      if ( .not. vmr_apply(ivmr) ) cycle
      call RF_VMR_Write( RF_VMR(region,ivmr), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end do

    ! lt output:
    if ( lt_apply ) then
      call RF_LT_Write( RF_LT(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if
    if ( lt2_apply ) then
      call RF_LT_Write( RF_LT2(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if

!     ! deposition velocities:
!     if ( depv_apply ) then
!       call RF_DEPV_Write( RF_DEPV(region), region, idate_f, status )
!       IF_NOTOK_RETURN(status=1)
!     end if

    status = 0

  END SUBROUTINE PDUMP_FILES_WRITE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PDUMP_FILES_WRITE2
  !
  ! !DESCRIPTION: write at end of time interval
  ! 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PDUMP_FILES_WRITE2( region, idate_f, status )
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    ! deposition fluxes:
!     if ( deps_apply ) then
!       call RF_DEPS_Write( RF_DEPS(region), region, idate_f, status )
!       IF_NOTOK_RETURN(status=1)
!     end if

    ! lt output:
    if ( lt_apply ) then
      call RF_LT_Write( RF_LT(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if
    if ( lt2_apply ) then
      call RF_LT_Write( RF_LT2(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! ok
    status = 0

  END SUBROUTINE PDUMP_FILES_WRITE2
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PDUMP_FILES_CLOSE
  !
  ! !DESCRIPTION: call done method of each output file.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PDUMP_FILES_CLOSE( region, status )
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)              ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)             ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !   31 Aug 2012 - Ph. Le Sager - switch closing order, since it was giving issues on some machine.
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer            ::  ivmr

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

!     if ( depv_apply ) then
!        call RF_DEPV_Done( RF_DEPV(region), status )
!        IF_NOTOK_RETURN(status=1)
!     end if

!     if ( deps_apply ) then
!        call RF_DEPS_Done( RF_DEPS(region), status )
!        IF_NOTOK_RETURN(status=1)
!     end if

    if ( lt2_apply ) then
       call RF_LT_Done( RF_LT2(region), region, status )
       IF_NOTOK_RETURN(status=1)
    end if
    
    if ( lt_apply ) then
      call RF_LT_Done( RF_LT(region), region, status )
      IF_NOTOK_RETURN(status=1)
    end if

    do ivmr = nvmr, 1, -1
       if ( .not. vmr_apply(ivmr) ) cycle
       call RF_VMR_Done( RF_VMR(region,ivmr), status )
       IF_NOTOK_RETURN(status=1)
    end do

    if ( tp_apply ) then
       call RF_TP_Done     ( RF_TP(region)     , status )
       IF_NOTOK_RETURN(status=1)
    end if

    if ( griddef_apply ) then
       call RF_GridDef_Done( RF_GridDef(region), status )
       IF_NOTOK_RETURN(status=1)
    end if

    status = 0

  end subroutine PDUMP_Files_Close
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_GRIDDEF_INIT
  !
  ! !DESCRIPTION:
  ! 
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE 1: Model horizontal grid definition 
  ! (longitude, latitude, size of gridbox [m2] ). 
  ! For documentation purposes, please also include the native vertical 
  ! grid definition from your model (hybrid level coefficients) and the 
  ! formula used to calculate pressure.
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine RF_GridDef_Init( RF, fdir, model, expid, region, status )
    !
    ! !USES:
    !
    use partools,  only : MPI_INFO_NULL, localComm
    use MeteoData, only : global_lli, levi
    !
    ! !OUTPUT PARAMETERS:
    !
    type(TPdumpFile_GridDef), intent(out)  ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)           ::  fdir
    character(len=*), intent(in)           ::  model
    character(len=*), intent(in)           ::  expid
    integer, intent(in)                    ::  region
    integer, intent(out)                   ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !   10 Jul 2012 - Ph. Le Sager - switch to MDF_NETCDF4
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    character(len=256)    ::  fname
    integer               ::  varid
    integer               ::  rtype

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

    call goLabel(rname)

    ! o open file

    ! write filename
    write (fname,'(a,"/",a,a,"_",a,"_",a,".nc")') &
         trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), 'griddef'

#ifdef MPI     
    ! overwrite existing files (clobber), provide MPI stuff:
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, &
                                  mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
    if (status/=0) then
       write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr
       write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr
       TRACEBACK; status=1; return
    end if
#else
    ! overwrite existing files (clobber)
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! o global attributes

    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title', 'model horizontal definition' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_author'     , trim(dataset_author)      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution'        , trim(institution)         , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version'    , trim(dataset_version)     , status)
    IF_NOTOK_MDF(fid=RF%ncid)


    ! o define dimensions

    call MDF_Def_Dim( RF%ncid, 'scalar', 1,                       RF%dimid_scalar , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lon',    global_lli(region)%nlon, RF%dimid_lon ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lat',    global_lli(region)%nlat, RF%dimid_lat ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lev',    levi%nlev,               RF%dimid_lev ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'levi',   levi%nlev+1,             RF%dimid_levi ,   status)
    IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Def_Dim( RF%ncid, 'time', NTS, RF%dimid_time , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Def_Dim( RF%ncid, 'datelen', 6, RF%dimid_datelen , status)
    !IF_NOTOK_MDF(fid=RF%ncid)

    ! o define variables
    rtype = MDF_FLOAT

    call MDF_Def_Var( RF%ncid, 'lon', rtype, (/RF%dimid_lon/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'longitude' ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'degrees_east' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lon = varid

    call MDF_Def_Var( RF%ncid, 'lat', MDF_FLOAT, (/RF%dimid_lat/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' ,      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'latitude' ,      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'degrees_north' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lat = varid

    !call MDF_Def_Var( RF%ncid, 'time', MDF_FLOAT, RF%dimid_time, varid , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    !IF_NOTOK_MDF(fid=RF%ncid)    
    !call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'time' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'units', 'days since 1950-01-01 00:00:00' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'calender', 'gregorian' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !RF%varid_time = varid
    
    !call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    !IF_NOTOK_MDF(fid=RF%ncid)    
    !call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'date' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !RF%varid_date = varid
    
    call MDF_Def_Var( RF%ncid, 'area', MDF_FLOAT, (/RF%dimid_lon,RF%dimid_lat/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'grid_cell_area' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'grid-cell area' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'm2' ,             status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_gridbox_area = varid

    call MDF_Def_Var( RF%ncid, 'a', MDF_FLOAT, (/RF%dimid_lev/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid sigma coordinate a coefficient' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a(k)*p0 + b(k)*ps(n,j,i)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'comment',       'bottom-up' ,                             status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_a = varid

    call MDF_Def_Var( RF%ncid, 'b', mdf_float, (/RF%dimid_lev/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid sigma coordinate a coefficient' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a(k)*p0 + b(k)*ps(n,j,i)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'comment',       'bottom-up' ,                             status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_b = varid

    call MDF_Def_Var( RF%ncid, 'a_bnds', mdf_float, (/RF%dimid_levi/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' ,                     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid sigma coordinate a coefficient for layer bounds' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                                      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_a_bnds = varid

    call MDF_Def_Var( RF%ncid, 'b_bnds', mdf_float, (/RF%dimid_levi/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' ,                     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid sigma coordinate a coefficient for layer bounds' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                                      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_b_bnds = varid

    call MDF_Def_Var( RF%ncid, 'p0', mdf_float, (/RF%dimid_scalar/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'reference pressure value' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_p0 = varid

    !status = pnf90_def_var( RF%ncid, 'ps', MDF_FLOAT, &
    !                         (/RF%dimid_lon,RF%dimid_lat,RF%dimid_time/), varid )
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'surface pressure' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'units', 'Pa' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !RF%varid_ps = varid

    !status = pnf90_def_var( RF%ncid, 'geo_height', MDF_FLOAT, &
    !              (/RF%dimid_lon,RF%dimid_lat,RF%dimid_lev,RF%dimid_time/), varid )
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'long_name', 'geopotential height' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'units', 'm' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !call MDF_Put_Att( RF%ncid, varid, 'comment', 'bottom-up; lower half level; top value implicit infinity' , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    !RF%varid_geo_height = varid

    !  o end defintion mode

    call MDF_EndDef( RF%ncid , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    ! no records written yet
    RF%trec = 0

    call goLabel() ; status = 0

  END SUBROUTINE RF_GRIDDEF_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  RF_GridDef_Write
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_GRIDDEF_WRITE( RF, region, status )
    !
    ! !USES:
    !
    use GO,        only : TDate, NewDate, rTotal, operator(-)
    use Grid,      only : AreaOper
    use MeteoData, only : global_lli, levi, sp_dat
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_GridDef), intent(inout)  ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                      ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                     ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !   10 Jul 2012 - Ph. Le Sager - switch to MDF_NETCDF4
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer             ::  imr, jmr, lmr
    real, allocatable   ::  ll(:,:)
    real                ::  time

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

    call goLabel(rname)

    ! grid size
    imr = global_lli(region)%nlon
    jmr = global_lli(region)%nlat
    lmr = levi%nlev

    ! next time record:
    RF%trec = RF%trec + 1

    ! o write data

    if ( RF%trec == 1 ) then

       ! lat/lon field:
       allocate( ll(imr,jmr) )

       call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg,  status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg,  status)
       IF_NOTOK_MDF(fid=RF%ncid)

       ll = 1.0
       call AreaOper( global_lli(region), ll, '*', 'm2', status )
       IF_NOTOK_RETURN(status=1)

       call MDF_Put_Var( RF%ncid, RF%varid_gridbox_area, ll , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_a, levi%fa , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_b, levi%fb , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_a_bnds, levi%a(0:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_b_bnds, levi%b(0:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_p0, (/1.0/) , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       deallocate( ll )

    end if
    
    !call MDF_Put_Var( RF%ncid, RF%varid_time, time, index=RF%trec , status)
    !IF_NOTOK_MDF(fid=RF%ncid)
    
    !call MDF_Put_Var( RF%ncid, RF%varid_date, reshape(real(idate_f),(/6,1/), status), &
    !                                       start=(/1,RF%trec/), count=(/6,1/) )
    !IF_NOTOK_MDF(fid=RF%ncid)
    
    !status = pnf90_put_var( RF%ncid, RF%varid_ps, &
    !                reshape(sp_dat(region)%data(1:imr,1:jmr,1:1),(/imr,jmr,1/)), &
    !                start=(/1,1,RF%trec/), count=(/imr,jmr,1/) )
    !IF_NOTOK_MDF(fid=RF%ncid)
    
    !status = pnf90_put_var( RF%ncid, RF%varid_geo_height, &
    !                 reshape(gph_dat(region)%data(1:imr,1:jmr,1:lmr),(/imr,jmr,lmr,1/)), &
    !                 start=(/1,1,1,RF%trec/), count=(/imr,jmr,lmr,1/) )
    !IF_NOTOK_MDF(fid=RF%ncid)

    
    call goLabel()
    status = 0

  END SUBROUTINE RF_GridDef_Write
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_GRIDDEF_DONE
  !
  ! !DESCRIPTION: close file-1
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_GridDef_Done( RF, status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_GridDef), intent(inout) :: RF
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,                  intent(out)   :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    call goLabel(rname)

    call MDF_Close( RF%ncid , status)
    IF_NOTOK_RETURN(status=1)

    call goLabel()

    status = 0

  END SUBROUTINE RF_GRIDDEF_DONE
  !EOC



  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE2: 3D field of monthly Model pressure [Pa] and temperature [K].
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_TP_INIT
  !
  ! !DESCRIPTION: file-2 : open and define var/att
  ! 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_TP_Init( RF, fdir, model, expid, region, idate_f, dhour, status )
    !
    ! !USES:
    !
    use partools,  only : MPI_INFO_NULL, localComm
    use MeteoData, only : global_lli, levi
    !
    ! !OUTPUT PARAMETERS:
    !
    type(TPdumpFile_TP), intent(out) ::  RF
    integer,             intent(out) ::  status
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)  ::  fdir
    character(len=*), intent(in)  ::  model
    character(len=*), intent(in)  ::  expid
    integer,          intent(in)  ::  region
    integer,          intent(in)  ::  idate_f(6)
    integer,          intent(in)  ::  dhour
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - switch to netcdf-4 thru MDF
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    character(len=256)    ::  fname
    integer               ::  varid, i1, i2, j1, j2

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

    call goLabel(rname)

    ! store arguments
    RF%dhour = dhour

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    n_tp_rec = GET_N_TIME_RECORDS( idate_f, dhour, mess='TP_Init' )

    if ( n_tp_rec == 0 ) then
       tp_apply = .false.
       status=0
       return
    end if

    
    ! o open file

    ! write filename
    write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') &
         trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), 'TP', idate_f(1:3)

    ! open, overwrite existing files (clobber)
#ifdef MPI     
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, &
                     mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
    if (status/=0) then
       write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr
       write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr
       TRACEBACK; status=1; return
    end if
#else
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! o global attributes

    call mdf_put_att( RF%ncid, MDF_GLOBAL, 'title', 'model pressure and temperature', status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, MDF_GLOBAL, 'dataset_author' , trim(dataset_author)  , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, MDF_GLOBAL, 'institution'    , trim(institution)     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, MDF_GLOBAL, 'dataset_version', trim(dataset_version) , status)
    IF_NOTOK_MDF(fid=RF%ncid)    

    ! o define dimensions
       
    call mdf_def_dim( RF%ncid, 'lon', global_lli(region)%nlon, RF%dimid_lon , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_def_dim( RF%ncid, 'lat', global_lli(region)%nlat, RF%dimid_lat , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_def_dim( RF%ncid, 'lev', levi%nlev,        RF%dimid_lev        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_def_dim( RF%ncid, 'time', n_tp_rec,        RF%dimid_time       , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_def_dim( RF%ncid, 'datelen', 6,            RF%dimid_datelen    , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    ! o define variables

    call mdf_def_var( RF%ncid, 'lon', MDF_FLOAT, (/RF%dimid_lon/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'longitude'       , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'longitude'       , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'degrees_east'    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lon = varid

    call mdf_def_var( RF%ncid, 'lat', MDF_FLOAT, (/RF%dimid_lat/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'latitude'        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'latitude'        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'degrees_north'   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lat = varid

    call mdf_def_var( RF%ncid, 'lev', MDF_FLOAT, (/RF%dimid_lev/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'level' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         '1'     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lev = varid

    call mdf_def_var( RF%ncid, 'time', MDF_FLOAT, (/RF%dimid_time/), varid              , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'time'                           , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'time'                           , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'days since 1950-01-01 00:00:00' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'calender',      'gregorian'                      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_time = varid

    allocate(RF%time(n_tp_rec))
    
    call mdf_def_var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'long_name', 'date and time' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_date = varid

    allocate(RF%date(6,n_tp_rec))
    
    call mdf_def_var( RF%ncid, 'ps', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'surface_air_pressure' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'surface pressure'     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'Pa'                   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_ps = varid

    call mdf_def_var( RF%ncid, 'orog', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'surface_altitude' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'surface altitude' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'm'                , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_orog = varid

    call mdf_def_var( RF%ncid, 'surface_temp', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'surface_temperature' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'surface temperature' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'K'                   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment', &
         '2m temperature from MARS archive or IFS model (grib 167, 2T)'      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_surface_temp = varid

    allocate( RF%data2d(i1:i2, j1:j2, n_tp_rec, 3) )


    call mdf_def_var( RF%ncid, 'geopotential', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'geopotential'    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'geopotential'    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'm2 s-2'          , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment',       'at mid levels'   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_geop = varid

    call mdf_def_var( RF%ncid, 'pressure', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'pressure'        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'pressure'        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'Pa'              , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment',       'at mid levels'   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_pressure = varid

    call mdf_def_var( RF%ncid, 'temp', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'air_temperature'        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'temperature'            , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'K'                      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment',       'bottom-up; full levels' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_temp = varid

    call mdf_def_var( RF%ncid, 'specific_humidity', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'specific_humidity' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'specific humidity' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'kg kg-1'           , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment', 'mass fraction of water vapor in moist air; (kg water)/(kg air)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_humid = varid

    call mdf_def_var( RF%ncid, 'U', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'eastward_wind'   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'zonal wind'      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'm s-1'           , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment', 'computed from mass fluxes through grid box boundaries' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_u = varid

    call mdf_def_var( RF%ncid, 'V', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'standard_name', 'northward_wind'  , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name',     'meridional wind' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',         'm s-1'           , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment', 'computed from mass fluxes through grid box boundaries' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_v = varid

    call mdf_def_var( RF%ncid, 'W', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'long_name', 'vertical wind velocity' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'units',     'm s-1'                  , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'comment', 'computed from mass fluxes through grid box boundaries' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_w = varid

    allocate( RF%data3d(i1:i2, j1:j2, levi%nlev, n_tp_rec, 7) )
    
    !  o end defintion mode

    call mdf_enddef( RF%ncid , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    !  o

    ! no records written yet
    RF%trec = 0

    call goLabel()

    ! ok
    status = 0

  END SUBROUTINE RF_TP_Init
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_TP_Write
  !
  ! !DESCRIPTION: store records, and if last time step write data to file
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_TP_Write( RF, region, idate_f, status )
    !
    ! !USES:
    !
    use Binas      , only : grav
    use Phys       , only : GeoPotentialHeight
    use Grid       , only : FPressure, HPressure
    use GO         , only : TDate, NewDate, rTotal, operator(-)
    use partools   , only : myid, root
    use MeteoData  , only : global_lli, lli, levi
    use MeteoData  , only : sp_dat, temper_dat, humid_dat, pu_dat, pv_dat, mfw_dat, gph_dat, oro_dat, t2m_dat
    use MeteoData  , only : m_dat
    use global_data, only : mass_dat
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_TP), intent(inout)     ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                    ::  region
    integer, intent(in)                    ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                   ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer           ::  i, j, l, i1, i2, j1, j2
    integer           ::  imr, jmr, lmr, klm
    real              ::  lev(levi%nlev)
    type(TDate)       ::  t, t0
    real              ::  time
    real, allocatable ::  field3d(:,:,:)
    real              ::  p_hlev(0:levi%nlev)

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

    ! for multiple of dhour only ...
    if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0)  ) then
      status=0; return
    end if

    call goLabel(rname)

    ! grid size
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=i2-i1+1
    jmr=j2-j1+1
    lmr = levi%nlev

    ! next time record:
    RF%trec = RF%trec + 1

    ! time since reftime:
    t0   = NewDate( time6=time_reftime6 )
    t    = NewDate( time6=idate_f )
    time = rTotal( t - t0, 'day' )
    
    if(okdebug)then
       write(gol,*) "RF_TP_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr
    end if
    
    ! o write data

    if ( RF%trec == 1 ) then

       call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       do l = 1, lmr
          lev(l) = real(l)
       end do
       call MDF_Put_Var( RF%ncid, RF%varid_lev, lev , status)
       IF_NOTOK_MDF(fid=RF%ncid)

    end if

    ! temporary storage for 3D fields
    allocate( field3d(i1:i2,j1:j2,1:lmr) ) ; field3d = 0.

    !-------- FILL DIAGNOSTIC ARRAYS
    RF%time(RF%trec)         = time
    RF%date(:,RF%trec)       = real(idate_f)
    RF%data2d(:,:,RF%trec,1) =  sp_dat(region)%data(i1:i2,j1:j2,1)
    RF%data2d(:,:,RF%trec,2) = oro_dat(region)%data(i1:i2,j1:j2,1)
    RF%data2d(:,:,RF%trec,3) = t2m_dat(region)%data(i1:i2,j1:j2,1)


    ! o  geopotential

    ! fill mid level geopotential:
    do j = j1, j2
       do i = i1, i2

          ! half level pressures
          call HPressure( levi, sp_dat(region)%data(i,j,1), p_hlev, status )
          IF_NOTOK_RETURN(status=1)

           ! mid level gph (m)
           call GeoPotentialHeight( lmr, p_hlev, temper_dat(region)%data(i,j,:), &
                humid_dat(region)%data(i,j,:), oro_dat(region)%data(i,j,1)/grav, &
                field3d(i,j,:) )    ! m
       end do
    end do

    ! multiply with gravity for correct unit:
    field3d = field3d * grav   ! m2/s2

    RF%data3d(:,:,:,RF%trec,1) =  field3d
    

    ! o  pressure

    ! fill mid level pressure
    call FPressure( levi, sp_dat(region)%data(i1:i2,j1:j2,1), field3d, status )
    IF_NOTOK_RETURN(status=1)

    RF%data3d(:,:,:,RF%trec,2) =  field3d

    ! o  temperature

    RF%data3d(:,:,:,RF%trec,3) =  temper_dat(region)%data(i1:i2,j1:j2,1:lmr)

    ! o  specific humidity

    RF%data3d(:,:,:,RF%trec,4) =  humid_dat(region)%data(i1:i2,j1:j2,1:lmr)

    ! o wind fields

    CALL UPDATE_HALO( dgrid(region), pu_dat(region)%data, pu_dat(region)%halo, status)
    IF_NOTOK_RETURN(status=1)

    CALL UPDATE_HALO( dgrid(region), pv_dat(region)%data, pv_dat(region)%halo, status)
    IF_NOTOK_RETURN(status=1)
    
    ! average U wind
    field3d = 0.5 * ( pu_dat(region)%data(i1-1:i2-1,j1:j2,1:lmr) + pu_dat(region)%data(i1:i2,j1:j2,1:lmr) ) &
         / m_dat(region)%data(i1:i2,j1:j2,1:lmr)    ! 1/s
    do j = j1, j2
       field3d(:,j,:) = field3d(:,j,:) * lli(region)%dx(j-j1+1)    ! m/s
    end do

    RF%data3d(:,:,:,RF%trec,5) =  field3d

    ! average V wind:
    field3d = 0.5 * ( pv_dat(region)%data(i1:i2,j1-1:j2-1,1:lmr) + pv_dat(region)%data(i1:i2,j1:j2,1:lmr) ) &
         / m_dat(region)%data(i1:i2,j1:j2,1:lmr)    ! 1/s
    field3d = field3d * lli(region)%dy   ! m/s

    RF%data3d(:,:,:,RF%trec,6) =  field3d

    ! from downward massflux to upward average W wind:
    field3d = 0.5 * ( mfw_dat(region)%data(i1:i2,j1:j2,0:lmr-1) + mfw_dat(region)%data(i1:i2,j1:j2,1:lmr) ) &
         / m_dat(region)%data(i1:i2,j1:j2,1:lmr)    ! 1/s
    do l = 1, lmr
       field3d(:,:,l) = - 1.0 * field3d(:,:,l) * &
            abs( gph_dat(region)%data(i1:i2,j1:j2,l+1) - gph_dat(region)%data(i1:i2,j1:j2,l) )   ! m/s
    end do

    RF%data3d(:,:,:,RF%trec,7) =  field3d


    
    !-------- WRITE ARRAYS
    if ( RF%trec == n_tp_rec ) then 

       ! time
       call MDF_Put_Var( RF%ncid, RF%varid_time, RF%time, status)!, start=(/1/), count=(/n_tp_rec/))
       IF_NOTOK_MDF(fid=RF%ncid)

       ! date
       call MDF_Put_Var( RF%ncid, RF%varid_date, RF%date, status )!, &
       !         start=(/1,1/), count=(/6,1/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! surface pressure
       call MDF_Put_Var( RF%ncid, RF%varid_ps, RF%data2d(:,:,:,1), status, start=(/i1,j1,1/), count=(/imr,jmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! orography (in m!)
       call MDF_Put_Var( RF%ncid, RF%varid_orog, RF%data2d(:,:,:,2), status, start=(/i1,j1,1/), count=(/imr,jmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! surface temperature = 2m temperature
       call MDF_Put_Var( RF%ncid, RF%varid_surface_temp, RF%data2d(:,:,:,3), status, start=(/i1,j1,1/) ) !, count=(/imr,jmr,1/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! geopotential
       call MDF_Put_Var( RF%ncid, RF%varid_geop, RF%data3d(:,:,:,:,1), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! pressure
       call MDF_Put_Var( RF%ncid, RF%varid_pressure, RF%data3d(:,:,:,:,2), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! temperature
       call MDF_Put_Var( RF%ncid, RF%varid_temp, RF%data3d(:,:,:,:,3), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! specific humidity
       call MDF_Put_Var( RF%ncid, RF%varid_humid, RF%data3d(:,:,:,:,4), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! winds
       call MDF_Put_Var( RF%ncid, RF%varid_u, RF%data3d(:,:,:,:,5), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Var( RF%ncid, RF%varid_v, RF%data3d(:,:,:,:,6), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Var( RF%ncid, RF%varid_w, RF%data3d(:,:,:,:,7), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,n_tp_rec/) )
       IF_NOTOK_MDF(fid=RF%ncid)
    end if
    
    ! Done
    deallocate( field3d )
    call goLabel()
    status = 0

  END SUBROUTINE RF_TP_Write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_TP_Done
  !
  ! !DESCRIPTION: close file #2
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine RF_TP_Done( RF, status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_TP),  intent(inout)  ::  RF
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,              intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    call goLabel(rname)

    call MDF_Close( RF%ncid , status)
    IF_NOTOK_RETURN(status=1)

    deallocate( rf%time, rf%date, rf%data2d, rf%data3d )

    call goLabel() ; status = 0

  end subroutine RF_TP_Done
  !EOC


  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE3: 3D fields for CO2 Volume Mixing Ratios
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_VMR_Init
  !
  ! !DESCRIPTION: open and define variables/attribute for file #3
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine RF_VMR_Init( RF, fdir, model, expid, filetype, region, &
                          idate_f, dhour,    tracer_names,  status  )
    !
    ! !USES:
    !
    use Binas,      only : xmair
    use GO,         only : goReadFromLine, goUpCase
    use chem_param, only : ntrace, names, ra
    use partools,   only : MPI_INFO_NULL, localComm
    use MeteoData,  only : global_lli, lli, levi, sp_dat
    use dims,       only : xbeg, xend, ybeg, yend, dx, dy, dz, xref, yref, zref
    use dims,       only : zbeg, zend
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_VMR), intent(inout) :: RF
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in) :: fdir
    character(len=*), intent(in) :: model
    character(len=*), intent(in) :: expid
    character(len=*), intent(in) :: filetype
    integer,          intent(in) :: region
    integer,          intent(in) :: idate_f(6)
    integer,          intent(in) :: dhour
    character(len=*), intent(in) :: tracer_names
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,              intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !   15 Apr 2014 - Ph. Le Sager - tropomi add-ons
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    character(len=256)    ::  fname, history, sysdate, model_meteo
    integer               ::  varid, i1, i2, j1, j2
    integer, dimension(8) ::  isysdate
    
    character(len=256)    ::  trnames
    character(len=8)      ::  trname, tmname
    integer               ::  k, itr, posend, pospoint
    integer               ::  imr, jmr, lmr, si, ei, ix, jy

    character(len=32)     ::  varname_spec
    character(len=5)      ::  zone

    character(len=64)     ::  cf_medium_stnd, cf_medium_long
    character(len=64)     ::  cf_enti_stnd, cf_enti_long, cf_enti_unit
    character(len=64)     ::  cf_spec_stnd, cf_spec_long
    character(len=4)      ::  cf_enti_type
    character(len=256)    ::  cf_name_stnd, cf_name_long, cf_name_unit
    character(len=512)    ::  comment
    character(len=6)      ::  csize

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

    call goLabel(rname)

    ! store arguments
    RF%dhour        = dhour
    RF%tracer_names = tracer_names

    ! size
    imr = global_lli(region)%nlon
    jmr = global_lli(region)%nlat
    lmr = levi%nlev

    ! number of time steps
    rf%n_rec = GET_N_TIME_RECORDS( idate_f, dhour, mess='VMR_Init' )

    ! degenerated cases (eg, very short runs)
    if ( rf%n_rec == 0 ) then
       rf%apply = .false.
       status=0
       return
    end if

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    

    ! set tracer index for requested tracers:
    write (gol,'("selected tracers for VMR output:")'); call goPr
    ! initialise RF
    RF%ntr  = 0
    RF%itr  = -1

    trnames = tracer_names
    do
      ! empty ?
      if ( len_trim(trnames) == 0 ) exit

      ! next number:
      if ( RF%ntr == ntrace ) then
        write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr
        TRACEBACK; status=1; return
      end if
      RF%ntr = RF%ntr + 1

      ! extract leading name:
      call goReadFromLine( trnames, trname, status, sep=' ' )
      IF_NOTOK_RETURN(status=1)

        ! convert to tm5 name:
        select case ( trim(strlowercase(trname)) )
        case default           ; tmname = trname
        end select

        ! --------------------------------
        ! NOy and M7 are special cases ...
        ! --------------------------------
        select case ( trim(strlowercase(tmname)) )

        case default
          ! --------------------------------
          ! `regular` constituents
          ! --------------------------------

          ! loop over all names:
          RF%itr(RF%ntr) = -1
          do itr = 1, ntrace
            ! case indendent match ?
            if ( goUpCase(trim(tmname)) == goUpCase(trim(names(itr))) ) then
              write (gol,'("  ",i3,"  ",a10,"  (",a10,") ",f12.4)') itr, trim(trname), trim(names(itr)), ra(itr); call goPr
              RF%itr(RF%ntr) = itr
              exit
            end if
          end do

        end select


        ! not found ?
        if ( RF%itr(RF%ntr) < 0 ) then
          write (gol,'("tracer name not supported:")'); call goPr
          write (gol,'("  list all     : ",a)') trim(tracer_names); call goPr
          write (gol,'("  list element : ",i3)') RF%ntr; call goPr
          write (gol,'("  pdump name   : ",a)') trim(trname); call goPr
          write (gol,'("  tm5 name     : ",a)') trim(tmname); call goPr
          write (gol,'("  tm5 tracers  : ")'); call goPr
          do itr = 1, ntrace
            write (gol,'("    ",i3," ",a)') itr, trim(names(itr)); call goPr
          end do
          TRACEBACK; status=1; return
        end if ! RF%itr 

      ! store pdump name:
      RF%name_tr(RF%ntr) = tmname

    end do



    ! empty file ?
    if ( RF%ntr < 1 ) then
      write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr
      TRACEBACK; status=1; return
    end if

    ! o open file

    ! write filename
    write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') &
         trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), trim(filetype), idate_f(1:3)

    ! open:
#ifdef MPI     
    ! overwrite existing files (clobber), provide MPI stuff:
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, &
                     mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
    if (status/=0) then
       write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr
       write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr
       TRACEBACK; status=1; return
    end if
#else
    ! overwrite existing files (clobber)
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! o global attributes

    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title'              , 'mixing ratios & concentrations' , status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution'        , trim(institution)                , status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version'    , trim(dataset_version)            , status )
    IF_NOTOK_MDF(fid=RF%ncid) 
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'file_version_number', trim(outfileversnr)              , status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'im'                 , imr                              , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                                         
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'jm'                 , jmr                              , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                                         
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'lm'                 , lmr                              , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                                         
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dx'                 , dx/xref(region)                  , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                                         
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dy'                 , dy/yref(region)                  , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                                         
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dz'                 , dz/zref(region)                  , status )
    IF_NOTOK_MDF(fid=RF%ncid)                 
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'xbeg'               , xbeg(region)                     , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'xend'               , xend(region)                     , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'ybeg'               , ybeg(region)                     , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'yend'               , yend(region)                     , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'zbeg'               , zbeg(region)                     , status )
    IF_NOTOK_MDF(fid=RF%ncid)                                
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'zend'               , zend(region)                     , status )
    IF_NOTOK_MDF(fid=RF%ncid)

    ! Meteo attribute
    if (trim(meteo_class)=='ei') then
       model_meteo='analysis (ERA-Interim)'
    elseif (trim(meteo_class)=='od') then
       model_meteo='forecast (IFS)'
    else
       write (gol,'("Meteo Model not known !")'); call goErr
       TRACEBACK; status=1; return
    endif
    
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'meteo_model', trim(model_meteo), status )
    IF_NOTOK_MDF(fid=RF%ncid)

    ! History attribute for audit trail: date, time of day, user name, program name
    call date_and_time(values=isysdate, zone=zone)

    write (sysdate, '(i4.4,"-",i2.2,"-",i2.2," ",i2.2,":",i2.2,":",i2.2," ",a)') &
         isysdate(1), isysdate(2), isysdate(3), isysdate(5), isysdate(6), isysdate(7), zone

    write(history,'("Created ",a," by ",a," with TM5.")') trim(sysdate),trim(dataset_author)

    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'history', trim(history), status )
    IF_NOTOK_MDF(fid=RF%ncid)     

       
    ! o define dimensions

    call MDF_Def_Dim( RF%ncid, 'lon',     imr,         RF%dimid_lon ,     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lat',     jmr,         RF%dimid_lat ,     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lev',     levi%nlev,   RF%dimid_lev ,     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'time',    rf%n_rec,    RF%dimid_time ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'levi',    levi%nlev+1, RF%dimid_levi ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'datelen', 6,           RF%dimid_datelen , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    ! o define variables

    call MDF_Def_Var( RF%ncid, 'lon', mdf_float, (/RF%dimid_lon/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude' ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'longitude' ,    status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'degrees_east' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lon = varid

    call MDF_Def_Var( RF%ncid, 'lat', mdf_float, (/RF%dimid_lat/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude' ,      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'latitude' ,      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'degrees_north' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lat = varid

    call MDF_Def_Var( RF%ncid, 'lev', mdf_float, (/RF%dimid_lev/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'level' ,                                       status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                           status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lev = varid

    call MDF_Def_Var( RF%ncid, 'time', mdf_double, (/RF%dimid_time/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time' ,                           status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'time' ,                           status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'days since 1950-01-01 00:00:00' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'calender',      'gregorian' ,                      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_time = varid

    allocate(RF%time(rf%n_rec))

    call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time'                          , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',     'year, month, day, hour, minute, second' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_date = varid

    allocate(RF%date(6,rf%n_rec))
    
    call MDF_Def_Var( RF%ncid, 'ps', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_time/), varid, status )
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'surface_air_pressure' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'surface pressure' ,     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'Pa' ,                   status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_ps = varid

    allocate( RF%sp(i1:i2, j1:j2, rf%n_rec) )

#ifdef tropomi
    call MDF_Def_Var( RF%ncid, 'temp', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'air_temperature'        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'temperature'            , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'K'                      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_put_att( RF%ncid, varid, 'comment',       'bottom-up; full levels' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_temp = varid	    

    allocate( RF%data3d_t(i1:i2, j1:j2, levi%nlev, rf%n_rec) )
#endif
    
    call MDF_Def_Var( RF%ncid, 'a_bnds', mdf_float, (/RF%dimid_levi/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' ,                     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid sigma coordinate a coefficient for layer bounds' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                                      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_a_bnds = varid

    call MDF_Def_Var( RF%ncid, 'b_bnds', mdf_float, (/RF%dimid_levi/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_coordinate' ,                     status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid sigma coordinate a coefficient for layer bounds' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '1' ,                                                      status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_b_bnds = varid


    ! loop over tracer to be written:
    do k = 1, RF%ntr

        ! ----------------------------
        ! setting defaults (gas phase)
        ! ----------------------------
        ! CF standard name for concentration/mixing ratio/column:
        cf_enti_stnd     = 'mole_fraction'
        cf_enti_unit     = 'mole mole-1'
        cf_enti_long     = 'volume mixing ratio'
        cf_medium_stnd   = 'in_air'
        cf_medium_long   = 'in humid air'

        RF%varid_type(k) = 'mixr'

        ! global tracer index
        itr = RF%itr(k)

        ! no comment yet
        comment = ''

        ! standard names from CF conventions:
        select case ( strlowercase(RF%name_tr(k)) )
        case ( 'co2'   )
          varname_spec = 'co2'
          cf_spec_stnd = 'carbon_dioxide'
          cf_spec_long = 'CO2'
        case default
          write (gol,'("do not know how to match tracer with CF standard names : ",a)') RF%name_tr(k); call goErr
          TRACEBACK; status=1; return
        end select


      ! define variable:
      call MDF_Def_Var( RF%ncid, trim(varname_spec), MDF_FLOAT, &
                        (/RF%dimid_lon,RF%dimid_lat,RF%dimid_lev,RF%dimid_time/), varid, status )
      IF_NOTOK_MDF(fid=RF%ncid)
      call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
      IF_NOTOK_MDF(fid=RF%ncid)    

      ! total names:
      cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)//'_'//trim(cf_medium_stnd)
      cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long)//' '//trim(cf_medium_long)
      cf_name_unit = trim(cf_enti_unit)


      ! write attributes:
      call MDF_Put_Att( RF%ncid, varid, 'standard_name', trim(cf_name_stnd) , status)
      IF_NOTOK_MDF(fid=RF%ncid)
      call MDF_Put_Att( RF%ncid, varid, 'long_name',     trim(cf_name_long) , status)
      IF_NOTOK_MDF(fid=RF%ncid)
      call MDF_Put_Att( RF%ncid, varid, 'units',         trim(cf_name_unit) , status)
      IF_NOTOK_MDF(fid=RF%ncid)

      ! moleweights; ra from chem_param is in g/mol .
      if ( itr <= ntrace .and. itr > 0 ) then
         call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', ra(itr)*1e3  , status)
         IF_NOTOK_MDF(fid=RF%ncid)
      else
         call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', -1.0         , status)
         IF_NOTOK_MDF(fid=RF%ncid)
      end if
      
      call MDF_Put_Att( RF%ncid   , varid, 'moleweight_air'   , xmair*1e3    , status)
      IF_NOTOK_MDF(fid=RF%ncid)

      call MDF_Put_Att( RF%ncid   , varid, 'moleweight_unit'  , 'kg mole-1'  , status)
      IF_NOTOK_MDF(fid=RF%ncid)

      if ( len_trim(comment) > 0 ) then
         call MDF_Put_Att( RF%ncid, varid, 'comment'          , trim(comment), status)
         IF_NOTOK_MDF(fid=RF%ncid)
      end if

      ! store varid
      RF%varid_tr(k) = varid

    end do

    ! storage
    allocate(rf%data3d(i1:i2,j1:j2,lmr,rf%n_rec,rf%ntr))

    !  o end defintion mode

    call MDF_EndDef( RF%ncid , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    !  o

    ! no records written yet
    RF%trec = 0

    call goLabel()
    status = 0

  END SUBROUTINE RF_VMR_Init
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_VMR_Write
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_VMR_Write( RF, region, idate_f, status )
    !
    ! !USES:
    !
    use Binas,       only : xmair
    use GO,          only : TDate, NewDate, rTotal, operator(-)
    use binas,       only : Rgas
    use chem_param,  only : ntrace, ntracet, fscale, ra
    use tracer_data, only : mass_dat, chem_dat
    use Grid,        only : FPressure
    use MeteoData,   only : global_lli, levi, m_dat, temper_dat, sp_dat
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_VMR), intent(inout)    ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                    ::  region
    integer, intent(in)                    ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                   ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !    2 Oct 2012 - Ph. Le Sager - adapted for lat-lon mpi decomp
    !                              - no more sub-regions available
    !
    ! !REMARKS:
    !  (1)
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer             ::  imr, jmr, lmr, i1, i2, j1, j2
    real, allocatable   ::  lev(:)
    integer             ::  l
    type(TDate)         ::  t, t0
    real                ::  time
    integer             ::  k, itr
    integer             ::  k_comp, itr_comp
    integer             ::  ims, ime, jms, jme, lms, lme
    integer             ::  gimr, gjmr, glmr

    real, allocatable   ::  compo_k(:,:,:)
    real, allocatable   ::  field_t(:,:,:)
    real, allocatable   ::  field_k(:,:,:)
    real, allocatable   ::  pres3d(:,:,:), pmx(:,:,:)
    integer             ::  numtrac
    integer             ::  listtrac(10)

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

    ! for multiple of dhour only ...
    if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0)  ) then
       status=0; return
    end if
   
    call goLabel(rname)

    ! grid sizes
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=i2-i1+1
    jmr=j2-j1+1
    lmr = levi%nlev
    
    gimr = global_lli(region)%nlon
    gjmr = global_lli(region)%nlat
    ! yet to change ?? 
    lms  = 1
    lme  = levi%nlev
    lmr  = levi%nlev
    glmr = levi%nlev

    ! next time record:
    RF%trec = RF%trec + 1
    if(okdebug)then
       write(gol,*) "RF_VMR_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr
    end if

    ! time since 1950-1-1 00:00
    t0   = NewDate( time6=time_reftime6 )
    t    = NewDate( time6=idate_f )
    time = rTotal( t - t0, 'day' )
    
    ! only once ...    
    if ( RF%trec == 1 ) then
       
       ! write longitudes:
       call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       ! write latitudes:
       call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       ! write level indices:
       allocate( lev(lmr) )
       do l = lms, lme
          lev(l) = real(l)
       end do
       call MDF_Put_Var( RF%ncid, RF%varid_lev, lev , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       deallocate(lev)

       ! As and Bs
       call MDF_Put_Var( RF%ncid, RF%varid_a_bnds, levi%a(0:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       call MDF_Put_Var( RF%ncid, RF%varid_b_bnds, levi%b(0:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)

    end if  ! first record


    RF%time(RF%trec)   = time
    RF%date(:,RF%trec) = real(idate_f)
    RF%sp(:,:,RF%trec) = sp_dat(region)%data(i1:i2,j1:j2,1)
#ifdef tropomi
    RF%data3d_t(:,:,:,RF%trec) = temper_dat(region)%data(i1:i2,j1:j2,1:lmr)
#endif    
    ! loop over all tracers to be written:
    do k = 1, RF%ntr

       ! global tracer index:
       itr = RF%itr(k)


          ! ---------
          ! transported or chemistry only ?
          ! ---------
          select case( itr )

          case( 1:ntracet ) 
                          
             ! ----------------------------------------------------
             ! distinguish between mixing ratios and concentrations
             ! ----------------------------------------------------
             select case( RF%varid_type(k) )

             case( 'conc' )
                ! write slab of concentrations
                !     m(trace)              pressure      xm(trace)
                ! C = -------- * fscale * -----------  *  ---------
                !      m(air)             temperature       Rgas 
!                 call MDF_Put_Var( RF%ncid, RF%varid_tr(k), &
!                                   reshape( mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr) / &
!                                            m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & 
!                                            pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & 
!                                            Rgas, (/imr,jmr,lmr,1/) ), &
!                                   status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) )

                rf%data3d(:,:,:, rf%trec, k) = mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr) / &
                                               m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & 
                                               pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & 
                                               Rgas
                
             case( 'mixr' )                
                ! write slab of volume mixing ratios
                !     m(trace)
                ! X = -------- * fscale 
                !      m(air)
!                 call MDF_Put_Var( RF%ncid, RF%varid_tr(k), &
!                                   reshape( mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr)/ &
!                                            m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr), &
!                                            (/imr,jmr,lmr,1/) ), &
!                                  status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) )

                rf%data3d(:,:,:, rf%trec, k) = mass_dat(region)%rm(i1:i2,j1:j2,lms:lme,itr)/ &
                                               m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr)
                
             case default 
                write (gol,'("no such unit type",a)') RF%varid_type(k); call goErr
                status=1
             end select

!             IF_NOTOK_MDF(fid=RF%ncid)
             

          ! ---------        
          case( ntracet+1:ntrace )
          ! ---------

             ! ----------------------------------------------------
             ! distinguish between mixing ratios and concentrations
             ! ----------------------------------------------------
             select case( RF%varid_type(k) )
             case( 'conc' )
                ! write slab of concentrations
                !     m(trace)              pressure      xm(trace)
                ! C = -------- * fscale * -----------  *  ---------
                !      m(air)             temperature       Rgas 
!                 call MDF_Put_Var( RF%ncid, RF%varid_tr(k), &
!                                   reshape( chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr) / &
!                                            m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & 
!                                            pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & 
!                                            Rgas, (/imr,jmr,lmr,1/) ), &
!                                            status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) )

                rf%data3d(:,:,:, rf%trec, k) = chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr) / &
                                               m_dat(region)%data(i1:i2,j1:j2,lms:lme) * xmair * 1.E-03 * & 
                                               pres3d(i1:i2,j1:j2,lms:lme) / temper_dat(region)%data(i1:i2,j1:j2,lms:lme) / & 
                                               Rgas
                
             case( 'mixr' )                
                ! write slab of volume mixing ratios
                !     m(trace)
                ! X = -------- * fscale 
                !      m(air)
!                 call MDF_Put_Var( RF%ncid, RF%varid_tr(k), &
!                                   reshape( chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr)/ &
!                                            m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr), &
!                                            (/imr,jmr,lmr,1/) ), &
!                                            status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) )

                rf%data3d(:,:,:, rf%trec, k) = chem_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr)/ &
                                               m_dat(region)%data(i1:i2,j1:j2,lms:lme) * fscale(itr)
                
             case default 
                write (gol,'("no such unit type",a)') RF%varid_type(k); call goErr
                status=1
             end select
                
             IF_NOTOK_MDF(fid=RF%ncid)

            ! -------------------
         case default 
            ! -------------------

            write (gol,'("strange tracer index requested : ",i6)') itr; call goErr
            TRACEBACK; status=1; return

         end select

   end do   ! tracer

   !----------------
   !    WRITE
   !----------------
   if ( RF%trec == rf%n_rec ) then

      call MDF_Put_Var( RF%ncid, RF%varid_time, rf%time, status)
      IF_NOTOK_MDF(fid=RF%ncid)

      call MDF_Put_Var( RF%ncid, RF%varid_date, rf%date, status)
      IF_NOTOK_MDF(fid=RF%ncid)

      ! surface presure
      call MDF_Put_Var( RF%ncid, RF%varid_ps, rf%sp, status, start=(/i1,j1,1/) )
      IF_NOTOK_MDF(fid=RF%ncid)

#ifdef tropomi      
      ! temperature
      call MDF_Put_Var( RF%ncid, RF%varid_temp, RF%data3d_t(:,:,:,:), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,RF%n_rec/) )
      IF_NOTOK_MDF(fid=RF%ncid)      
#endif
      ! vmr
      do k = 1, RF%ntr

         call MDF_Put_Var( RF%ncid, RF%varid_tr(k), RF%data3d(:,:,:,:,k), status, start=(/i1,j1,1,1/) )         
         IF_NOTOK_MDF(fid=RF%ncid)

      end do

   end if


   !----------------
   ! DONE
   !----------------
   call goLabel()
   status = 0

  END SUBROUTINE RF_VMR_Write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_VMR_Done
  !
  ! !DESCRIPTION: close file #3
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_VMR_Done( RF, status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_VMR), intent(inout)  ::  RF
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,              intent(out)    ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    call goLabel(rname)

    call MDF_Close( RF%ncid, status )
    IF_NOTOK_RETURN(status=1)

    deallocate(rf%date, rf%time, rf%sp, rf%data3d )

#ifdef tropomi
    deallocate(rf%data3d_t)
#endif
    
    call goLabel() ; status = 0

  END SUBROUTINE RF_VMR_Done
  !EOC


  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE: 2D LT output
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_LT_Init
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine RF_LT_Init( RF, fdir, model, expid, filetype, region, &
                         idate_f, local_time, tracer_names, status )
    !
    ! !USES:
    !
    use Binas,      only : xmair
    use GO,         only : goReadFromLine, goUpCase
    use GO,         only : NewDate
    use dims,       only : im, jm
    use chem_param, only : ntrace, names, ra
    use partools,   only : MPI_INFO_NULL, localComm
    use MeteoData,  only : global_lli, levi, sp_dat, Set
    !
    ! !OUTPUT PARAMETERS:
    !
    type(TPdumpFile_LT), intent(out)       ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*),  intent(in)   ::  fdir
    character(len=*),  intent(in)   ::  model
    character(len=*),  intent(in)   ::  expid
    character(len=*),  intent(in)   ::  filetype
    integer,           intent(in)   ::  region
    integer,           intent(in)   ::  idate_f(6)
    integer,           intent(in)   ::  local_time
    character(len=*),  intent(in)   ::  tracer_names
    integer,           intent(out)  ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    character(len=256)    ::  fname
    integer               ::  varid
    integer               ::  imr, jmr, lmr

    character(len=256)    ::  trnames
    character(len=8)      ::  trname, tmname
    character(len=3)      ::  cwavel
    integer               ::  k, itr, i1, i2, j1, j2

    character(len=32)     ::  varname, varname_enti, varname_spec

    character(len=64)     ::  cf_medium_stnd, cf_medium_long
    character(len=64)     ::  cf_enti_stnd, cf_enti_long, cf_enti_unit
    character(len=64)     ::  cf_spec_stnd, cf_spec_long
    character(len=256)    ::  cf_name_stnd, cf_name_long, cf_name_unit
    character(len=512)    ::  comment

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

    call goLabel(rname)

    ! store arguments
    RF%local_time       = local_time
    RF%tracer_names     = tracer_names

    ! set tracer index for requested tracers:
    write (gol,'("selected tracers for LT output:")'); call goPr
    RF%ntr = 0
    RF%itr = -1

    trnames = tracer_names
    do
       ! empty ?
       if ( len_trim(trnames) == 0 ) exit
       
       ! next number:
       if ( RF%ntr == ntrace ) then
          write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr
          TRACEBACK; status=1; return
       end if
       RF%ntr = RF%ntr + 1

       ! extract leading name:
       call goReadFromLine( trnames, trname, status, sep=' ' )
       IF_NOTOK_RETURN(status=1)

          ! convert to tm5 name:
          select case ( trim(strlowercase(trname)) )
          case default           ; tmname = trname
          end select

          ! NOy is a special ...
          select case ( trim(strlowercase(tmname)) )
          case default
             ! loop over all names:
             RF%itr(RF%ntr) = -1
             do itr = 1, ntrace
                ! case indendent match ?
                if ( goUpCase(trim(tmname)) == goUpCase(trim(names(itr))) ) then
                   write (gol,'("  ",i3,"  ",a10,"  (",a10,") ",f12.4)') itr, trim(trname), trim(names(itr)), ra(itr); call goPr
                   RF%itr(RF%ntr) = itr
                   exit
                end if
             end do
          end select      ! not found ?


          if ( RF%itr(RF%ntr) < 0 ) then
             write (gol,'("tracer name not supported:")'); call goPr
             write (gol,'("  list all     : ",a)') trim(tracer_names); call goPr
             write (gol,'("  list element : ",i3)') RF%ntr; call goPr
             write (gol,'("  pdump name   : ",a)') trim(trname); call goPr
             write (gol,'("  tm5 name     : ",a)') trim(tmname); call goPr
             write (gol,'("  tm5 tracers  : ")'); call goPr
             do itr = 1, ntrace
                write (gol,'("    ",i3," ",a)') itr, trim(names(itr)); call goPr
             end do
             TRACEBACK; status=1; return
          end if

       ! store pdump name:
       RF%name_tr(RF%ntr) = trname

    end do
    
    ! empty file ?
    if ( RF%ntr < 1 ) then
       write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr
       TRACEBACK; status=1; return
    end if

    ! grid size
    imr = global_lli(region)%nlon
    jmr = global_lli(region)%nlat
    lmr = levi%nlev


    ! o open file

    ! write filename
    write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".nc")') &
         trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), trim(filetype), idate_f(1:3)

    ! open:
#ifdef MPI     
    ! overwrite existing files (clobber), provide MPI stuff:
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status, &
                     mpi_comm=localComm, mpi_info=MPI_INFO_NULL )
    if (status/=0) then
       write (gol,'("from creating NetCDF4 file for writing in parallel;")'); call goErr
       write (gol,'("MDF module not compiled with netcdf4_par support ?")'); call goErr
       TRACEBACK; status=1; return
    end if
#else
    ! overwrite existing files (clobber)
    call MDF_Create( trim(fname), MDF_NETCDF4, MDF_REPLACE, RF%ncid, status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    ! o global attributes

    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'title'              , 'local time output'      , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_author'     , trim(dataset_author)     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution'        , trim(institution)        , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_version'    , trim(dataset_version)    , status)
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'file_version_number', trim(outfileversnr)      , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    ! o define dimensions

    call MDF_Def_Dim( RF%ncid, 'lon'    , global_lli(region)%nlon, RF%dimid_lon    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lat'    , global_lli(region)%nlat, RF%dimid_lat    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'lev'    , levi%nlev       , RF%dimid_lev    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'time'   , 1               , RF%dimid_time   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Def_Dim( RF%ncid, 'datelen', 6               , RF%dimid_datelen, status)
    IF_NOTOK_MDF(fid=RF%ncid)

    ! o define variables

    call MDF_Def_Var( RF%ncid, 'lon', mdf_float, (/RF%dimid_lon/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'longitude'   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name'    , 'longitude'   , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units'        , 'degrees_east', status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lon = varid

    call MDF_Def_Var( RF%ncid, 'lat', mdf_float, (/RF%dimid_lat/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'latitude'     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name'    , 'latitude'     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units'        , 'degrees_north', status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lat = varid

    call MDF_Def_Var( RF%ncid, 'lev', mdf_float, (/RF%dimid_lev/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'atmosphere_hybrid_sigma_pressure_coordinate'     , status)
    IF_NOTOK_MDF(fid=RF%ncid)                                                                            
    call MDF_Put_Att( RF%ncid, varid, 'long_name'    , 'level'                                           , status)
    IF_NOTOK_MDF(fid=RF%ncid)                                                                            
    call MDF_Put_Att( RF%ncid, varid, 'units'        , '1'                                               , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'p(n,k,j,i) = a_bnds(k)*p0 + b_bnds(k)*ps(n,j,i)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_lev = varid

    call MDF_Def_Var( RF%ncid, 'time', mdf_float, (/RF%dimid_time/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'time'                          , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name'    , 'time'                          , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units'        , 'days since 1950-01-01 00:00:00', status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'calender'     , 'gregorian'                     , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_time = varid

    call MDF_Def_Var( RF%ncid, 'date', MDF_FLOAT, (/RF%dimid_datelen,RF%dimid_time/), varid , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'date and time' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units', 'year, month, day, hour, minute, second' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_date = varid


    call MDF_Def_Var( RF%ncid, 'ps', MDF_FLOAT, &
         (/RF%dimid_lon,RF%dimid_lat,RF%dimid_time/), varid, status )
    call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
    IF_NOTOK_MDF(fid=RF%ncid)    
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'standard_name', 'surface_air_pressure', status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name'    , 'surface pressure'    , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units'        , 'Pa'                  , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_ps = varid

    ! CF standard name for medium:
    cf_medium_stnd = 'in_air' ; cf_medium_long = 'in humid air'

    ! loop over tracer to be written:
    do k = 1, RF%ntr

          ! global tracer index
          itr = RF%itr(k)

          ! ~~ local time species info

          ! CF standard name for concentration/mixing ratio/column:
          cf_enti_stnd = 'mole_fraction'
          cf_enti_unit = 'mole mole-1'
          cf_enti_long = 'volume mixing ratio'

          ! start of dataset name:
          varname_enti = 'dry'

          ! no comment yet
          comment = ''

          ! standard names from CF conventions:
          select case ( RF%name_tr(k) )
          case ( 'CO2', 'co2'   )
             varname_spec = 'co2'
             cf_spec_stnd = 'carbon_dioxide'
             cf_spec_long = 'CO2'
          case default
             write (gol,'("do not know how to match tracer with CF standard names : ",a)') RF%name_tr(k); call goPr
             TRACEBACK; status=1; return
          end select

       ! define variable:
       call MDF_Def_Var( RF%ncid, trim(varname_spec), MDF_FLOAT, &
                        (/RF%dimid_lon,RF%dimid_lat,RF%dimid_lev,RF%dimid_time/), varid, status )
       IF_NOTOK_MDF(fid=RF%ncid)    

       call MDF_Var_Par_Access( RF%ncid, varid, access_mode, status )
       IF_NOTOK_MDF(fid=RF%ncid)

       ! total names:
       cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)//'_'//trim(cf_medium_stnd)
       cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_long)//' '//trim(cf_medium_long)
       cf_name_unit = trim(cf_enti_unit)      

       ! write attributes:
       call MDF_Put_Att( RF%ncid,    varid, 'standard_name',     trim(cf_name_stnd) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Att( RF%ncid,    varid, 'long_name',         trim(cf_name_long) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Att( RF%ncid,    varid, 'units',             trim(cf_name_unit) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       
       if ( itr <= ntrace .and. itr > 0 ) then
          call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', ra(itr)*1e3 ,        status)
          IF_NOTOK_MDF(fid=RF%ncid)
       else
          call MDF_Put_Att( RF%ncid, varid, 'moleweight_tracer', -1.0 ,               status)
          IF_NOTOK_MDF(fid=RF%ncid)
       end if

       call MDF_Put_Att( RF%ncid,    varid, 'moleweight_air',    xmair*1e3 ,          status)
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Att( RF%ncid,    varid, 'moleweight_unit',   'kg mole-1' ,        status)
       IF_NOTOK_MDF(fid=RF%ncid)

       if ( len_trim(comment) > 0 ) then
          call MDF_Put_Att( RF%ncid, varid, 'comment',           trim(comment) ,      status)
          IF_NOTOK_MDF(fid=RF%ncid)
       end if
       
       ! store varid
       RF%varid_tr(k) = varid
       
    end do

    !  o end defintion mode

    call MDF_EndDef( RF%ncid , status)
    IF_NOTOK_MDF(fid=RF%ncid)

    ! no records written yet
    RF%trec = 0


    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    
    allocate(RF%accu   (i1:i2, j1:j2, 1:lmr, RF%ntr)) ; RF%accu = 0
    allocate(RF%naccu  (i1:i2, RF%ntr ))              ; RF%naccu = 0
    allocate(RF%p_accu (i1:i2, j1:j2  ))              ; RF%p_accu = 0
    allocate(RF%np_accu(i1:i2         ))              ; RF%np_accu = 0

    call goLabel()

    status = 0

  END SUBROUTINE RF_LT_Init
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_LT_Write
  !
  ! !DESCRIPTION: does not write anything, just get
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_LT_Write( RF, region, idate_f, status )
    !
    ! !USES:
    !
    use GO,          only : TDate, NewDate, Set, iTotal, rTotal, operator(-), wrtgol
    use chem_param,  only : ntrace, ntracet, fscale
    use tracer_data, only : mass_dat, chem_dat
    use MeteoData,   only : global_lli, levi, m_dat, sp_dat
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_LT),   intent(inout)   ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                    ::  region
    integer, intent(in)                    ::  idate_f(6)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                   ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer             ::  imr, jmr, lmr, gimr, i1, i2, j1, j2
    real, allocatable   ::  lev(:)
    real, allocatable   ::  field_out(:,:,:)
    real, allocatable   ::  field_out_b(:,:)
    integer             ::  l, ls, le
    type(TDate)         ::  t, t0
    real                ::  time
    real                ::  dt_sec
    integer             ::  i, j, k, itr, itau, loctim, gridboxtimestep
    integer             ::  iloctim,itautoday,ilon

    integer             ::  icomp, itr_loc, ncells, window

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

    ! for multiple of dhour only ...
    !    if ( (modulo(idate_f(4),RF%dhour)/=0) .or. any(idate_f(5:6)/=0)  ) then
    !      status=0; return
    !    end if

    call goLabel(rname)

    ! grid size
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=i2-i1+1
    jmr=j2-j1+1
    gimr = global_lli(region)%nlon
!    gjmr = global_lli(region)%nlat
    lmr = levi%nlev


    ! next time record:
    RF%trec = RF%trec + 1
    if(okdebug)then
       write(gol,*) "RF_LT_Write - idate_f(6), RF%trec=", idate_f, RF%trec; call goPr
    end if

    ! grid index offsets for GMT and local time
    loctim=RF%local_time
    if( loctim < 0 ) loctim=loctim+24*3600

    ! time since 1950-1-1 00:00
    t0 = NewDate( time6=time_reftime6 )
    t  = NewDate( time6=idate_f )
    call SET( t, hour=0, min=0, sec=0 )
    time = rTotal( t - t0, 'day' ) + loctim / 86400.

    !
    ! ~~ time, grid
    !
    
    ! only once ...    
    if ( RF%trec == 1 ) then

       ! write longitudes:
       call MDF_Put_Var( RF%ncid, RF%varid_lon, global_lli(region)%lon_deg , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       ! write latitudes:
       call MDF_Put_Var( RF%ncid, RF%varid_lat, global_lli(region)%lat_deg , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       ! write level indices:
       allocate( lev(lmr) )
       do l = 1, lmr
          lev(l) = real(l)
       end do
       call MDF_Put_Var( RF%ncid, RF%varid_lev, lev , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       deallocate(lev)

       ! time:
       call MDF_Put_Var( RF%ncid, RF%varid_time, (/time/) , status, start=(/RF%trec/))
       IF_NOTOK_MDF(fid=RF%ncid)

       ! date:
       call MDF_Put_Var( RF%ncid, RF%varid_date, reshape(real(idate_f),(/6,1/)), status, &
            start=(/1,1/), count=(/6,1/) )
       IF_NOTOK_MDF(fid=RF%ncid)

    end if  ! first record
    
    ! 
    ! local time 
    !
    if ( RF%trec > 1 ) then ! do not accumulate fields on 00:00

       ! grid index offsets for GMT and local time
       loctim=RF%local_time
       if( loctim < 0 ) loctim=loctim+24*3600


       gridboxtimestep=24*3600/gimr
       itau = idate_f(4)*3600+idate_f(5)*60+idate_f(6)


       itautoday= nint(real(mod(itau,24*3600)*gimr)/real(24*3600))
       iloctim  = nint(real(loctim           *gimr)/real(24*3600))

       ! determine longitude index wrt Greenwich from difference (local time - GMT)
       ! also process neigboring longitudes (i-2 , i-1 , i , i+1 , i+2) depending on 
       ! number of longitudinal grid cells
       ncells = ceiling( gimr / 24. )
       window = ceiling( ncells / 2. ) 
       do ilon = 1, ncells
          i = 1 + mod( gimr + gimr/2 + iloctim - itautoday + (ilon - window),gimr )
          if (i .ge. i1 .and. i.le. i2) then

             RF%p_accu(i,j1:j2)= RF%p_accu(i,j1:j2)+sp_dat(region)%data(i,j1:j2,1)
             RF%np_accu(i)= RF%np_accu(i)+1


             ! loop over tracers to be written:
             do k = 1, RF%ntr

                ! global tracer index:
                itr = RF%itr(k)


                ! transported or chemistry only ?
                if ( (itr >= 1) .and. (itr <= ntracet) ) then

                   RF%accu(i,j1:j2,1:lmr,k)= RF%accu(i,j1:j2,1:lmr,k)+&
                        (mass_dat(region)%rm(i,j1:j2,1:lmr,itr)/ &
                        m_dat(region)%data(i,j1:j2,1:lmr))*fscale(itr)

                   RF%naccu(i,k)=RF%naccu(i,k)+1

                else if ( (itr >= ntracet+1) .and. (itr <= ntrace) ) then

                   RF%accu(i,j1:j2,1:lmr,k)= RF%accu(i,j1:j2,1:lmr,k)+&
                        (chem_dat(region)%rm(i,j1:j2,1:lmr,itr)/ &
                        m_dat(region)%data(i,j1:j2,1:lmr))*fscale(itr)

                   RF%naccu(i,k)=RF%naccu(i,k)+1

                end if


             enddo

          endif
       enddo

    endif ! do not accumulate fields on 00:00


    call goLabel(); status = 0

  END SUBROUTINE RF_LT_Write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_LT_Done
  !
  ! !DESCRIPTION: write final data, then close file #4
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_LT_Done( RF, region, status )
    !
    ! !USES:
    !
    use MeteoData, only : global_lli, levi
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_LT), intent(inout) ::  RF
    !
    ! !INPUT PARAMETERS:
    !
    integer,             intent(in)    ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,             intent(out)   ::  status
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retro -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !                              - move averaging & writing here
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len =*), parameter ::  rname = mname//'/RF_LT_Done'
    integer                      ::  imr, jmr
    real, allocatable            ::  field_out(:,:,:)
    real, allocatable            ::  field_out_b(:,:)
    integer                      ::  i, ls, le, k, itr, i1, i2, j1, j2, lmr

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

    call goLabel(rname)

    !---------------------
    ! average & write data
    !---------------------    
    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=i2-i1+1
    jmr=j2-j1+1
    lmr = levi%nlev

    allocate(field_out_b(i1:i2,j1:j2)); field_out_b = 0.0

    do i = i1, i2
       if (RF%np_accu(i).gt.0) then
          field_out_b(i,:) =RF%p_accu(i,:)/RF%np_accu(i)
       endif
    enddo

    call MDF_Put_Var( RF%ncid, RF%varid_ps, reshape(field_out_b(i1:i2,j1:j2), &
                      (/imr,jmr,1/) ), status, start=(/i1,j1,1/), count=(/imr,jmr,1/) )
    IF_NOTOK_MDF(fid=RF%ncid)
    
    deallocate(field_out_b)


    TRACERS: do k = 1, RF%ntr

       ! global tracer index:
       itr = RF%itr(k)

       if ( (itr >= 1) .and. (itr <= ntrace) ) then

          ! normalize fields, if necessary
          allocate(field_out(i1:i2,j1:j2,1:lmr)); field_out = 0.0
          do i = i1,i2
             if (RF%naccu(i,k).gt.0) then
                field_out(i,:,1:lmr) =RF%accu(i,:,1:lmr,k)/RF%naccu(i,k)
             endif
          enddo

          ! write fields:
          call MDF_Put_Var( RF%ncid, RF%varid_tr(k)              , &
                            reshape(field_out(i1:i2,j1:j2,1:lmr) , &
                            (/imr,jmr,lmr,1/) )                  , &
                            status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,1/) )
          IF_NOTOK_MDF(fid=RF%ncid)
          deallocate(field_out)

       endif

    end do TRACERS


    !---------------------    
    ! DONE
    !---------------------
    call MDF_Close( RF%ncid , status)
    IF_NOTOK_RETURN(status=1)

    deallocate(RF%accu)
    deallocate(RF%naccu)
    deallocate(RF%p_accu)
    deallocate(RF%np_accu)

    call goLabel() ; status = 0

  END SUBROUTINE RF_LT_Done
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    strlowercase
  !
  ! !DESCRIPTION: 
  !
  !   This function returns a copy of the input string 'struppercase' with all
  !   letters changed to lowercase. All other characters remain unchanged.
  !\\
  !\\
  ! !INTERFACE:
  !
  FUNCTION strlowercase(struppercase)
    !
    ! !USES:
    !
    IMPLICIT NONE
    !
    ! !INPUT PARAMETERS:
    !
    CHARACTER(LEN=*), INTENT(IN)     :: struppercase
    !
    ! !RETURN VALUE:
    !
    CHARACTER(LEN=LEN(struppercase)) :: strlowercase
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(LEN=1)                 :: u
    INTEGER                          :: i,j
    strlowercase = struppercase
    DO i=1,LEN(struppercase)
      u = struppercase(i:i)
      j = IACHAR(u)
      IF(j < 65 .OR. j > 90) CYCLE
      strlowercase(i:i) = ACHAR(j+32)
    END DO
    !-------------------------------------------------------------------------------
  END FUNCTION STRLOWERCASE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  struppercase
  !
  ! !DESCRIPTION: 
  !
  !   This function returns a copy of the input string 'struppercase' with all
  !   letters changed to lowercase. All other characters remain unchanged.
  !\\
  !\\
  ! !INTERFACE:
  !
  FUNCTION STRUPPERCASE(strlowercase)
    !
    ! !USES:
    !
    IMPLICIT NONE
    !
    ! !INPUT PARAMETERS:
    !
    CHARACTER(LEN=*), INTENT(IN)     :: strlowercase
    !
    ! !RETURN VALUE:
    !
    CHARACTER(LEN=LEN(strlowercase)) :: struppercase
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk -
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    CHARACTER(LEN=1)                 :: u
    INTEGER                          :: i,j
    struppercase = strlowercase
    DO i=1,LEN(strlowercase)
      u = strlowercase(i:i)
      j = IACHAR(u)
      IF(j < 97 .OR. j > 122) CYCLE
      struppercase(i:i) = ACHAR(j-32)
    END DO
    !-------------------------------------------------------------------------------
  END FUNCTION STRUPPERCASE
  !EOC

END MODULE USER_OUTPUT_PDUMP