ecearth3/sources/tm5mp/proj/cb05/user_output_pdump.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
!
! If the macro (cpp) "tropomi" is used, then the temperature and extra attributes added to the vmr (tracers) datasets.
!
! 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
! tropomi only:
!    output.pdump.tropomi.tm5version     :  v4
!    output.pdump.tropomi.institution    :  KNMI
!    output.pdump.tropomi.tm5reference   :  Huijnen et al., ACP
!    output.pdump.tropomi.authoremail    :  Doe@john.com
!    output.pdump.tropomi.datasetname    :  "S5P_AUX_CTMFCT" or "S5P_AUX_CTMANA"
!
!    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          :  vmr3
!    output.pdump.vmr.001.dhour          :  3
!    output.pdump.vmr.001.tracers        :  SO2 NOy CH4 OH HNO3 PAN H2O2 Radon Lead
!
!    output.pdump.vmr.002.apply          :  T
!    output.pdump.vmr.002.fname          :  vmr1
!    output.pdump.vmr.002.dhour          :  1
!    output.pdump.vmr.002.tracers        :  O3 O3s CO NO2 NO CH2O
!
!    output.pdump.vmr.003.apply          :  F
!    output.pdump.vmr.003.fname          :  vmra
!    output.pdump.vmr.003.dhour          :  3
!    output.pdump.vmr.003.tracers        : SO4 NO3_A BC BCS POM SS1_N SS1_M SS2_N SS2_M SS3_N SS3_M DUST2_N DUST2_M DUST3_N DUST3_M
!
!    output.pdump.lt.apply               :  T
!    output.pdump.lt.tracers             :  O3
!    output.pdump.lt.localtime           :  2
!    
!    output.pdump.lt2.apply              :  F
!    output.pdump.lt2.tracers            :  
!    output.pdump.lt2.localtime          :  
! 
!    output.pdump.depositions.apply      :  F
!    output.pdump.depositions.dhour      :  3
!    output.pdump.depositions.tracers    :  O3 HNO3 NO NO2 H2O2 CH2O PAN CO NH3 NH4 SO2 NOy
!
!    output.pdump.depvels.apply          :  F
!    output.pdump.depvels.dhour          :  3
!    output.pdump.depvels.tracers        :  O3 HNO3 NO NO2 H2O2 CH2O PAN CO NH3 NH4 SO2
!
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_PDUMP
  !
  ! !USES:
  !
  use partools,     only : isRoot
  use GO,           only : gol, goPr, goErr, goLabel
  use GO,           only : TDate, IncrDate, NewDate
  use GO,           only : operator(+), SystemDate, Get
  use dims,         only : nregions, idatee, idatei, okdebug, nread
  use chem_param,   only : ntrace
  use chem_param,   only : iNOx, iHNO3, iPAN, iOrgNtr
#ifdef with_m7
  use chem_param,   only : iNO3_a
  use chem_param,   only : iSO4nus, iSO4ais, iSO4acs, iSO4cos
  use chem_param,   only : iBCais,  iBCacs,  iBCcos,  iBCaii
  use chem_param,   only : iPOMais, iPOMacs, iPOMcos, iPOMaii
  use chem_param,   only : iDUacs,  iDUcos,  iDUaci,  iDUcoi
  use chem_param,   only : iSScos,  iSSacs
#endif
  
  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' 
  !
  ! NOy is not a standard tracer field, but sum of some transported tracers:  
  !   NOx HNO3 PAN orgntr NO3_a
  ! where NOx is the sum of short lived tracers:
  !   NOx = NO + NO2 + NO3 + HNO4 + 2*N2O5 
  !
#ifdef with_m7
  integer, parameter  ::  iNOy         = ntrace + 1    
  integer, parameter  ::  nNOyt        = 5
  integer, parameter  ::  iNOyt(nNOyt) = (/ iNOx, iHNO3, iNO3_a, iPAN, iOrgNtr /)

  integer, parameter  ::  iSO4         = ntrace + 2
  integer, parameter  ::  nSO4t        = 4
  integer, parameter  ::  iSO4t(nSO4t) = (/ iSO4nus, iSO4ais, iSO4acs, iSO4cos /)

  integer, parameter  ::  iBC          = ntrace + 3
  integer, parameter  ::  nBCt         = 4
  integer, parameter  ::  iBCt(nBCt)   = (/ iBCais, iBCacs, iBCcos, iBCaii /)

  integer, parameter  ::  iPOM         = ntrace + 4
  integer, parameter  ::  nPOMt        = 4
  integer, parameter  ::  iPOMt(nPOMt) = (/ iPOMais, iPOMacs, iPOMcos, iPOMaii /)

  integer, parameter  ::  iSS          = ntrace + 5
  integer, parameter  ::  nSSt         = 2
  integer, parameter  ::  iSSt(nSSt)   = (/ iSSacs, iSScos /)

  integer, parameter  ::  iDU          = ntrace + 6
  integer, parameter  ::  nDUt         = 4
  integer, parameter  ::  iDUt(nDUt)   = (/ iDUacs, iDUcos, iDUaci, iDUcoi /)
#else
  integer, parameter  ::  iNOy         = ntrace + 1    
  integer, parameter  ::  nNOyt        = 4
  integer, parameter  ::  iNOyt(nNOyt) = (/ iNOx, iHNO3, iPAN, iOrgNtr /)
#endif 
  !
  ! !PRIVATE TYPES:
  !
  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
     real                ::  dhour
     integer             ::  dsec
    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_ps
    integer             ::  varid_a_bnds, varid_b_bnds
    integer             ::  ntr
    integer             ::  itr(ntrace)
    character(len=8)    ::  name_tr(ntrace)
#ifdef with_m7
    logical             ::  lpmx(ntrace)
    real                ::  sizepmx(ntrace)
#endif
    integer             ::  varid_tr(ntrace)
    character(len=4)    ::  varid_type(ntrace)
    real, allocatable   ::  data3d(:,:,:,:,:)
    real, allocatable   ::  sp(:,:,:)
    real, allocatable   ::  time(:)
    real, allocatable   ::  date(:,:)
    real, allocatable   ::  data3d_t(:,:,:,:)
    integer             ::  varid_temp    
#ifdef tropomi
    integer             ::  varid_hyai, varid_hybi, varid_hyam, varid_hybm
    integer             ::  varid_hgt
    integer             ::  varid_ltropo
    real, allocatable   ::  data2d_hgt(:,:)
    integer, allocatable::  data2d_ltropo(:,:,:)
#endif
  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(:)
#ifdef with_m7
    logical             ::  laod(ntrace)
    real                ::  wavel(ntrace)
#endif
  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
#ifdef tropomi
  character(len=256)              ::  tropomi_authoremail, tropomi_tm5_reference, tropomi_institution
  character(len=256)              ::  tropomi_tm5_version, tropomi_dataset_name
  character(len=15)               ::  tropomi_date_start, tropomi_date_stop, tropomi_date_create
#endif

  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(:)
  real, 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,
  !                                  handle aerosol tracers and M7 
  !   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
  !   8 October 2014 - H. Eskes - changes in tropomi output (based on the "tropomi" macro) 
  !
  ! !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.
  !   (5) Switch in nstep for DEPS data should work for full days. Not tested
  !         for partial days.
  !
  ! !TODO:
  !   - test with M7 tracers. Which ones?
  !   - in LT_WRITE    : AOD if m7 needs to be coded
  !   - in RF_VMR_INIT : match tracer with CF standard names for some aerosols
  !                      (dust,...)
  !
  !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, dsec, 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 (timestart basically)
    integer, intent(in)  :: dsec     ! time step for timeseries in sec (should divide 24*3600, be divided by ndyn/2)
    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"
    !   15 Jul 2014 - Ph. Le Sager - works with seconds instead of hours
    !
    ! !REMARKS:
    !  - dynamic timestep cannot be LARGER than timestep of timeseries, with notable exception
    !      of dynamic timestep = 2*timeseries_timestep. 
    !
    ! !TODO:
    !  - 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)*3600+date(5)*60+date(6)  ! 0, unless start of the run is not at 00:00:00
    if (deps) delta=delta + nread          ! one DYNAMIC time step done to output something

    if (modulo(delta,dsec)==0) then
       is=delta/dsec
    else
       is=(delta+dsec)/dsec
    end if
    
    ! End index for daily file (nread=dynamic time step read from rc)
    ie = (24*3600 - nread/2)  / dsec

    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*3600/dsec
    end if

    ! Case of "last day stopping before midnite". (Need testing for DEPS)
    if (lastday) ie=(idatee(4)*3600+idatee(5)*60+idatee(6)-nread/2)/dsec
    
    ! 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
       write(gol,*) "GET_N_TIME_RECORDS - date, dsec, nread ", date, dsec, nread                 ; call goPr
       write(gol,*) "GET_N_TIME_RECORDS - idateE ", idatee                                       ; 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, dsec_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)         ::  dsec_min  ! smallest timeseries period in sec
    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.
    
    ! lowest time frequency in sec
    dsec_min = 999999
    
    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)

#ifdef tropomi
    call ReadRc( rcF, 'output.pdump.tropomi.tm5version', tropomi_tm5_version , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.tropomi.institution', tropomi_institution , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.tropomi.tm5reference', tropomi_tm5_reference , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.tropomi.authoremail', tropomi_authoremail , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.pdump.tropomi.datasetname', tropomi_dataset_name , status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! 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

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

    ! 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. )
          call Set(    gph_dat(region), status, used=.true. )   ! used to compute surface altitude
       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)

             ! here is the catch: fractional hour for step should divide 3600 
             dsec_min = min( dsec_min, int(vmr_dhour(ivmr)*3600) )
             
             ! 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
#ifdef with_budgets
      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)
      dsec_min = min( dsec_min, deps_dhour*3600)
#else
      write(gol,*) "timeseries of deposition fluxes requires using 'with_budget' macro" ; call goErr
      status=1 ; TRACEBACK ; return
#endif
    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
#ifdef with_budgets
      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)
      dsec_min = min( dsec_min, depv_dhour*3600)
#else
      write(gol,*) "timeseries of deposition velocities requires using 'with_budget' macro" ; call goErr
      status=1 ; TRACEBACK ; return
#endif
    end if

    ! no files open yet
    curr_day = -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 )
    !
    ! !USES:
    !
    use dims,                 only : itaur
    use datetime,             only : tau2date
    !
    ! !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
    integer,dimension(6)        :: idate_f

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

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

      ! write end of interval DEPS data (requires that DEPS nstep is calculated with .false. -see RF_DEPS_Init-)
      call tau2date(itaur(region),idate_f)
      call PDUMP_Files_Write2( region, idate_f, status )
      IF_NOTOK_RETURN(status=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

#ifdef with_budgets
    ! 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
#endif

    ! 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

#ifdef with_budgets
    ! deposition velocities:
    if ( depv_apply ) then
      call RF_DEPV_Write( RF_DEPV(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if
#endif
    
    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 -------------------------------

#ifdef with_budgets
    ! deposition fluxes:
    if ( deps_apply ) then
      call RF_DEPS_Write( RF_DEPS(region), region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end if
#endif
    
    ! 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 -------------------------------

#ifdef with_budgets
    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
#endif
    
    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*3600, 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 O3, CO, CH4, ...  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 : PAR_BROADCAST, 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)
    real,             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
    !    8 Oct 2014 - H. Eskes     - 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

    integer, dimension(6) ::  idate_f_end, idate_create
    type(TDate)           ::  date_f_end, date_create
    
    ! --- begin -------------------------------------

    call goLabel(rname)

    ! store arguments
    RF%dhour        = dhour
    RF%dsec         = int(dhour*3600.)
    RF%tracer_names = tracer_names
    
    ! Test that dsec is multiple of dynamic-step/2 (nread in sec)
    if (((RF%dsec*2)/nread < 1).or.(modulo(RF%dsec,nread/2)/=0))then
       write(gol,*) "timeseries timestep should be a multiple of (dynamic_timestep)/2"; call goErr
       TRACEBACK; status=1; return
    end if
    
    ! 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, rf%dsec, 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
#ifdef with_m7
    RF%lpmx    = .false.
    RF%sizepmx = -1.0
#endif
    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)

#ifdef with_m7


      ! ---------------------------
      ! check for PMx 
      ! ---------------------------
      if( strlowercase(trname(1:2)) == 'pm' ) then 

        RF%lpmx(RF%ntr) = .true. 
        RF%itr (RF%ntr) = -1
        ! paste size to real 
        read(trname(3:len_trim(trname)), * )  RF%sizepmx(RF%ntr)

      else 

#endif         

        ! convert to tm5 name:
        select case ( trim(strlowercase(trname)) )
        case ( 'hcho'        ) ; tmname = 'CH2O'
        case ( 'rn', 'radon' ) ; tmname = 'Rn222'
        case ( 'pb', 'lead'  ) ; tmname = 'Pb210'
        case default           ; tmname = trname
        end select

        ! --------------------------------
        ! NOy and M7 are special cases ...
        ! --------------------------------
        select case ( trim(strlowercase(tmname)) )
        case( 'noy' )
          ! defined as ntrace+1
          RF%itr(RF%ntr) = iNOy
          write (gol,'("    *  ",a10)') trim(trname); call goPr
#ifdef with_m7
        case( 'tso4' ) 
          ! defined as ntrace+2
          RF%itr(RF%ntr) = iSO4
          write (gol,'("    *  ",a10)') trim(trname); call goPr
        case( 'tbc' ) 
          ! defined as ntrace+3
          RF%itr(RF%ntr) = iBC
          write (gol,'("    *  ",a10)') trim(trname); call goPr
        case( 'tpom' ) 
          ! defined as ntrace+4
          RF%itr(RF%ntr) = iPOM
          write (gol,'("    *  ",a10)') trim(trname); call goPr
        case( 'tss' ) 
          ! defined as ntrace+5
          RF%itr(RF%ntr) = iSS
          write (gol,'("    *  ",a10)') trim(trname); call goPr
        case( 'tdu' ) 
          ! defined as ntrace+6
          RF%itr(RF%ntr) = iDU
          write (gol,'("    *  ",a10)') trim(trname); call goPr

#endif

        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 

#ifdef with_m7
      end if ! pmx
#endif

      ! 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

#ifdef tropomi
    ! define start/stop of output, and run date
    date_f_end = NewDate( time6=idate_f ) + IncrDate(hour=24)
    call Get( date_f_end, time6=idate_f_end )
    if ( isRoot ) then
       date_create = SystemDate()
       call Get( date_create, time6=idate_create )
    endif
    call PAR_BROADCAST(idate_create, status)
    IF_NOTOK_RETURN(status=1)

    date_create = SystemDate()
    call Get( date_create, time6=idate_create ) 
    write (tropomi_date_start, '(i4.4,i2.2,i2.2,"T",i2.2,i2.2,i2.2)') idate_f
    write (tropomi_date_stop,  '(i4.4,i2.2,i2.2,"T",i2.2,i2.2,i2.2)') idate_f_end
    write (tropomi_date_create,'(i4.4,i2.2,i2.2,"T",i2.2,i2.2,i2.2)') idate_create

    ! write filename according to tropomi convention
    write (fname,'(a,"/",a,"_",a,"_",a,".nc")') &
         trim(fdir), trim(tropomi_dataset_name), tropomi_date_start, tropomi_date_stop
#else
    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)
#endif

    ! 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

#ifdef tropomi
    ! H. Eskes: Extra attributes for TROPOMI 
    ! Conventions = "CF-1.6"
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'Conventions', 'CF-1.6' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! validity_start = "20132305T120000" (zoals in filenaam)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'validity_start', tropomi_date_start , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! validity_stop = "20132405T000000"
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'validity_stop', tropomi_date_stop , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! creation_date = "20142909T124905"
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'creation_date', tropomi_date_create , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! version = TM5 version string.
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'version', trim(tropomi_tm5_version) , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! institution = "KNMI"
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'institution', trim(tropomi_institution) , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! reference = TM5 reference (journal article or so)
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'reference', trim(tropomi_tm5_reference) , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! contact = email address of volunteer.
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'contact', trim(tropomi_authoremail) , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    ! dataset_name = "S5P_NRTI_AUX_CTMFCT" of "S5P_OFFL_AUX_CTMANA"
    call MDF_Put_Att( RF%ncid, MDF_GLOBAL, 'dataset_name', trim(tropomi_dataset_name) , status)
    IF_NOTOK_MDF(fid=RF%ncid)
#endif

    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)=='ea') then
       model_meteo='reanalysis (ERA5)'
    elseif (trim(meteo_class)=='od') then
       model_meteo='forecast (IFS)'
    elseif (trim(meteo_class)=='ifs10') then
       model_meteo='EC-Earth (ifs 10L)'
     elseif (trim(meteo_class)=='ifs34') then
       model_meteo='EC-Earth (ifs 34L)'
     elseif (trim(meteo_class)=='ifs62') then
       model_meteo='EC-Earth (ifs 62L)'
    elseif (trim(meteo_class)=='ifs91') then
       model_meteo='EC-Earth (ifs 91L)'
    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, 'levi',    levi%nlev+1, RF%dimid_levi ,    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, '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
    
#ifdef tropomi
    call MDF_Def_Var( RF%ncid, 'hyai', 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, 'units',         'Pa' ,                                          status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid A coefficient at layer interfaces' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_hyai = varid
#else
    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
#endif

#ifdef tropomi
    call MDF_Def_Var( RF%ncid, 'hybi', 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, 'units',         '1' ,                                           status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid B coefficient at layer interfaces' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_hybi = varid
#else
    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
#endif
    
#ifdef tropomi
    call MDF_Def_Var( RF%ncid, 'hyam', 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, 'units',         'Pa' ,                                          status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid A coefficient at layer midpoints' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_hyam = varid
    call MDF_Def_Var( RF%ncid, 'hybm', 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, 'units',         '1' ,                                          status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name', 'hybrid B coefficient at layer midpoints' ,        status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_hybm = varid
#endif

    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)
#ifdef tropomi
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'hybrid level at layer midpoints' ,             status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         'level' ,                                       status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call mdf_put_att( RF%ncid, varid, 'positive',      'down'  , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula', 'hyam hybm (mlev=hyam+hybm*ps)' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'formula_terms', 'ap: hyam b: hybm ps: ps' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
#else
    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)
#endif
    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 )
    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

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

#ifndef tropomi
    call MDF_Def_Var( RF%ncid, 't', 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

#ifdef tropomi

    ! Extra temperature field output
    ! with compression - crash
    !call MDF_Def_Var( RF%ncid, 't', MDF_FLOAT, (/RF%dimid_lon, RF%dimid_lat, RF%dimid_lev, RF%dimid_time/), varid, status, compression=1, deflate_level=4)
    call MDF_Def_Var( RF%ncid, 't', 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) )

    ! Extra surface elevation output, retrieved from GPH (meteo.f90) and g0 (binas.f90) following WGS84?

    call MDF_Def_Var( RF%ncid, 'surface_altitude', 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', 'surface_altitude'             , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'surface altitude of TM5 grid' , 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',       'ECMWF interpolated orography' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_hgt = varid

    allocate( RF%data2d_hgt(i1:i2, j1:j2) )

    ! Extra tropopause level output, retrieved from GPH and temperature (meteo.f90)

    call MDF_Def_Var( RF%ncid, 'tropopause_layer_index', MDF_INT, (/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', 'tropopause_layer_index'             , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'long_name',     'index of the highest model layer in the troposphere' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid, 'units',         '-'                            , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_put_att( RF%ncid, varid, 'comment',       'Based on WMO temperature gradient method' , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_ltropo = varid

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

#endif

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

#ifdef with_m7
      if( RF%lpmx(k) ) then 

        ! get diameter 
        write(csize,'(F5.1)') RF%sizepmx(k) 

        ! remove leading blanks
        csize    = adjustl(csize) 
        pospoint = index(csize,'.')
        posend   = len_trim(csize)

        ! CF standard name for concentration/mixing ratio/column:
        RF%varid_type(k) = 'conc'
        varname_spec     = 'pm'//csize(1:pospoint-1)//'p'//csize(pospoint+1:posend)
        cf_spec_stnd     = 'particulate_matter_'//trim(csize)
        cf_spec_long     = 'particulate matter diameter le '//trim(csize)//' micrometers'
        cf_enti_stnd     = 'concentration'
        cf_enti_unit     = 'kg m-3 '
        cf_enti_long     = 'mass per volume'

      else 
#endif       

        ! ----------------------------
        ! setting defaults (gas phase)
        ! ----------------------------
        ! CF standard name for concentration/mixing ratio/column:
        cf_enti_stnd     = 'mole_fraction'
#ifdef tropomi
        cf_enti_unit     = '1'
#else
        cf_enti_unit     = 'mole mole-1'
#endif
        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 ( 'co'   )
          varname_spec = 'co'
          cf_spec_stnd = 'carbon_monoxide'
          cf_spec_long = 'CO'
        case ( 'o3'   )
          varname_spec = 'o3'
          cf_spec_stnd = 'ozone'
          cf_spec_long = 'O3'
        case ( 'o3s'   )
          varname_spec = 'o3s'
          cf_spec_stnd = 'ozone_from_stratosphere'
          cf_spec_long = 'O3s'
        case ( 'no'  )
          varname_spec = 'no'
          cf_spec_stnd = 'nitrogen_monoxide'
          cf_spec_long = 'NO'
        case ( 'no2'  )
          varname_spec = 'no2'
          cf_spec_stnd = 'nitrogen_dioxide'
          cf_spec_long = 'NO2'
        case ( 'noy'  )
          varname_spec = 'noy'
          cf_spec_stnd = 'nitrogen_oxides'
          cf_spec_long = 'NOy'
          comment      = 'NOy = NOx + HNO3 + PAN + org.ntr., '// &
               'with NOx = NO + NO2 + NO3 + HNO4 + N2O5'
        case ( 'ch2o', 'choh' )
          varname_spec = 'ch2o'
          cf_spec_stnd = 'formaldehyde'
          cf_spec_long = 'CH2O'
        case ( 'so2'  )
          varname_spec = 'so2'
          cf_spec_stnd = 'sulfur_dioxide'
          cf_spec_long = 'SO2'
        case( 'h2so4' )
          varname_spec = 'h2so4'
          cf_spec_stnd = 'sulfuric_acid_g'
          cf_spec_long = 'H2SO4 (g)'
!!$        case ( 'so4'  )
!!$          varname_spec = 'so4'
!!$          cf_spec_stnd = 'sulfate_as_sulfate_dry_aerosol'
!!$          cf_spec_long = 'SO4'
        case ( 'ch4'  )
          varname_spec = 'ch4'
          cf_spec_stnd = 'methane'
          cf_spec_long = 'CH4'
        case ( 'oh'  )
          varname_spec = 'oh'
          cf_spec_stnd = 'hydroxyl_radical'
          cf_spec_long = 'OH'
        case ( 'h2o2'  )
          varname_spec = 'h2o2'
          cf_spec_stnd = 'hydrogen_peroxide'
          cf_spec_long = 'H2O2'
        case ( 'hno3'  )
          varname_spec = 'hno3'
          cf_spec_stnd = 'nitric_acid'
          cf_spec_long = 'HNO3'

        case ( 'hno4'  )
          varname_spec = 'hno4'
          cf_spec_stnd = 'peroxonitric_acid'
          cf_spec_long = 'HNO4'
        case ( 'n2o5'  )
          varname_spec = 'n2o5'
          cf_spec_stnd = 'nitrogen_pentoxide'
          cf_spec_long = 'N2O5'
        case ( 'par'  )
          varname_spec = 'par'
          cf_spec_stnd = 'paraffinic_carbon_atoms'
          cf_spec_long = 'PAR'
        case ( 'eth'  )
          varname_spec = 'eth'
          cf_spec_stnd = 'ethylene'
          cf_spec_long = 'ETH'
        case ( 'ole'  )
          varname_spec = 'ole'
          cf_spec_stnd = 'olefinic_carbon_bonds'
          cf_spec_long = 'OLE'
        case ( 'ald2'  )
          varname_spec = 'ald2'
          cf_spec_stnd = 'acetaldehyde_and_higher_aldehydes'
          cf_spec_long = 'ALD2'
        case ( 'mgly'  )
          varname_spec = 'mgly'
          cf_spec_stnd = 'methylglyoxal'
          cf_spec_long = 'MGLY'
        case ( 'isop'  )
          varname_spec = 'isop'
          cf_spec_stnd = 'isoprene'
          cf_spec_long = 'ISOP'
        case ( 'nh3'  )
          varname_spec = 'nh3'
          cf_spec_stnd = 'ammonia'
          cf_spec_long = 'NH3'
        case ( 'ORGNTR','orgntr'  )
          varname_spec = 'orgntr'
          cf_spec_stnd = 'organic_nitrate'
          cf_spec_long = 'ORGNTR'
        case ( 'pan'  )
          varname_spec = 'pan'
          cf_spec_stnd = 'peroxyacetyl_nitrate'
          cf_spec_long = 'PAN'

        case ( 'terp'  )
          varname_spec = 'terp'
          cf_spec_stnd = 'terpene'
          cf_spec_long = 'TERP'
        case ( 'elvoc'  )
          varname_spec = 'elvoc'
          cf_spec_stnd = 'extremely low volatile OC'
          cf_spec_long = 'ELVOC'
        case ( 'svoc'  )
          varname_spec = 'svoc'
          cf_spec_stnd = 'semi volatile OC'
          cf_spec_long = 'SVOC'

        case ( 'rn', 'radon', 'rn222'  )
          varname_spec = 'rn'
          cf_spec_stnd = 'radon'
          cf_spec_long = 'Rn'
        case ( 'pb', 'lead', 'pb210'  )
          varname_spec = 'pb'
          cf_spec_stnd = 'lead'
          cf_spec_long = 'Pb'

#ifdef with_m7
          ! Sulphate
        case( 'tso4' )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'so4'
          cf_spec_stnd     = 'total_sulphate_aerosol'
          cf_spec_long     = 'SO4'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
          ! Black Carbon
        case( 'tbc' )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'bc'
          cf_spec_stnd     = 'total_black_carbon_aerosol'
          cf_spec_long     = 'BC'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
          ! Particulate Organic Matter
        case( 'tpom' )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'pom'
          cf_spec_stnd     = 'total_particulate_organic_matter_aerosol'
          cf_spec_long     = 'POM'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
          ! Sea Salt
        case( 'tss' )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'ss'
          cf_spec_stnd     = 'total_sea_salt_aerosol'
          cf_spec_long     = 'SS'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
          ! Dust
        case( 'tdu' )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'du'
          cf_spec_stnd     = 'total_dust_aerosol'
          cf_spec_long     = 'SS'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
          ! Nucleation Soluble (nus): number, SO4
        case ( 'nus_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'nus_n'
          cf_spec_stnd     = 'number_wet_nucleation'
          cf_spec_long     = 'Number_nus'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'so4nus'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'so4nus'
          cf_spec_stnd     = 'sulphate_wet_nucleation'
          cf_spec_long     = 'SO4_nus'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'soanus'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'soanus'
          cf_spec_stnd     = 'SOA_wet_nucleation'
          cf_spec_long     = 'SOA_nus'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
          ! Aitken Soluble (ais): number, SO4, BC, POM
        case ( 'ais_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'ais_n'
          cf_spec_stnd     = 'number_wet_aitken'
          cf_spec_long     = 'Number_ais'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'so4ais'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'so4ais'
          cf_spec_stnd     = 'sulphate_wet_aitken'
          cf_spec_long     = 'SO4_ais'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'bcais'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'bcais'
          cf_spec_stnd     = 'black_carbon_wet_aitken'
          cf_spec_long     = 'BC_ais'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'pomais'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'pomais'
          cf_spec_stnd     = 'particulate_organic_matter_wet_aitken'
          cf_spec_long     = 'POM_ais'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'soaais'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'soaais'
          cf_spec_stnd     = 'SOA_dry_Aitken'
          cf_spec_long     = 'SOA_ais'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'

          ! Accumulation Soluble (acs): number, SO4, BC, POM, SS, DU
        case ( 'acs_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'acs_n'
          cf_spec_stnd     = 'number_wet_accumulation'
          cf_spec_long     = 'Number_acs'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'so4acs'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'so4acs'
          cf_spec_stnd     = 'sulphate_wet_accumulation'
          cf_spec_long     = 'SO4_acs'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'bcacs'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'bcacs'
          cf_spec_stnd     = 'black_carbon_wet_accumulation'
          cf_spec_long     = 'BC_acs'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'pomacs'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'pomacs'
          cf_spec_stnd     = 'particulate_organic_matter_wet_accumulation'
          cf_spec_long     = 'POM_acs'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'ssacs'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'ssacs'
          cf_spec_stnd     = 'seasalt_wet_accumulation'
          cf_spec_long     = 'SS_acs'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'duacs'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'duacs'
          cf_spec_stnd     = 'dust_wet_accumulation'
          cf_spec_long     = 'DU_acs'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'soaacs'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'soaacs'
          cf_spec_stnd     = 'SOA_dry_Accumulation'
          cf_spec_long     = 'SOA_acs'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
 
          ! Coarse Soluble (cos): number, SO4, BC, POM, SS, DU
        case ( 'cos_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'cos_n'
          cf_spec_stnd     = 'number_wet_coarse'
          cf_spec_long     = 'Number_cos'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'so4cos'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'so4cos'
          cf_spec_stnd     = 'sulphate_wet_coarse'
          cf_spec_long     = 'SO4_cos'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'bccos'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'bccos'
          cf_spec_stnd     = 'black_carbon_wet_coarse'
          cf_spec_long     = 'BC_cos'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'pomcos'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'pomcos'
          cf_spec_stnd     = 'particulate_organic_matter_wet_coarse'
          cf_spec_long     = 'POM_cos'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'sscos'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'sscos'
          cf_spec_stnd     = 'seasalt_wet_coarse'
          cf_spec_long     = 'SS_cos'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'ducos'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'ducos'
          cf_spec_stnd     = 'dust_wet_coarse'
          cf_spec_long     = 'DU_cos'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'soacos'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'soacos'
          cf_spec_stnd     = 'SOA_dry_coarse'
          cf_spec_long     = 'SOA_cos'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'

          ! Aitken Insoluble (aii): number, BC, POM
        case ( 'aii_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'aii_n'
          cf_spec_stnd     = 'number_dry_aitken'
          cf_spec_long     = 'Number_aii'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'bcaii'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'bcaii'
          cf_spec_stnd     = 'black_carbon_dry_aitken'
          cf_spec_long     = 'BC_aii'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'pomaii'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'pomaii'
          cf_spec_stnd     = 'particulate_organic_matter_dry_aitken'
          cf_spec_long     = 'POM_aii'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'soaaii'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'soaaii'
          cf_spec_stnd     = 'SOA_dry_Aitken'
          cf_spec_long     = 'SOA_aii'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'

          ! Accumulation Insoluble (aci): number, DU
        case ( 'aci_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'aci_n'
          cf_spec_stnd     = 'number_dry_accumulation'
          cf_spec_long     = 'Number_aci'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'duaci'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'duaci'
          cf_spec_stnd     = 'dust_dry_accumulation'
          cf_spec_long     = 'DU_aci'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'

          ! Coarse Insoluble (coi): number, DU
        case ( 'coi_n'  )
          RF%varid_type(k) = 'numb'
          varname_spec     = 'coi_n'
          cf_spec_stnd     = 'number_dry_coarse'
          cf_spec_long     = 'Number_coi'
          cf_enti_stnd     = 'number'
          cf_enti_unit     = '1.'
          cf_enti_long     = ''
        case ( 'ducoi'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'ducoi'
          cf_spec_stnd     = 'dust_dry_coarse'
          cf_spec_long     = 'DU_coi'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
#endif

        case ( 'nh4'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'nh4'
          cf_spec_stnd     = 'ammonium_as_ammonium_dry_aerosol'
          cf_spec_long     = 'NH4'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
        case ( 'no3_a'  )
          RF%varid_type(k) = 'conc'
          varname_spec     = 'no3'
          cf_spec_stnd     = 'nitrate_as_nitrate_dry_aerosol'
          cf_spec_long     = 'NO3'
          cf_enti_stnd     = 'concentration'
          cf_enti_unit     = 'kg m-3 '
          cf_enti_long     = 'mass per volume'
!!$        case ( 'bc'  )
!!$          varname_spec = 'bc'
!!$          cf_spec_stnd = 'black_carbon_dry_aerosol'
!!$          cf_spec_long = 'BC'
!!$        case ( 'BCS', 'bcs'  )
!!$          varname_spec = 'bcs'
!!$          cf_spec_stnd = 'hydrophilic_black_carbon_dry_aerosol'
!!$          cf_spec_long = 'BC(aq)'
!!$        case ( 'POM', 'pom'  )
!!$          varname_spec = 'om'
!!$          cf_spec_stnd = 'organic_carbon_as_particulate_organic_matter_dry_aerosol'
!!$          cf_spec_long = 'OM'
!!$        case ( 'SS1_N', 'ss1_n' )
!!$          varname_spec = 'ss1_n'
!!$          cf_spec_stnd = 'seasalt_dry_aerosol_mode1_number'
!!$          cf_spec_long = 'SS1_n'
!!$        case ( 'SS1_M', 'ss1_m' )
!!$          varname_spec = 'ss1_m'
!!$          cf_spec_stnd = 'seasalt_dry_aerosol_mode1_mass'
!!$          cf_spec_long = 'SS1_m'
!!$        case ( 'SS2_N', 'ss2_n' )
!!$          varname_spec = 'ss2_n'
!!$          cf_spec_stnd = 'seasalt_dry_aerosol_mode2_number'
!!$          cf_spec_long = 'SS2_n'
!!$        case ( 'SS2_M', 'ss2_m' )
!!$          varname_spec = 'ss2_m'
!!$          cf_spec_stnd = 'seasalt_dry_aerosol_mode2_mass'
!!$          cf_spec_long = 'SS2_m'
!!$        case ( 'SS3_N', 'ss3_n' )
!!$          varname_spec = 'ss3_n'
!!$          cf_spec_stnd = 'seasalt_dry_aerosol_mode3_number'
!!$          cf_spec_long = 'SS3_n'
!!$        case ( 'SS3_M', 'ss3_m' )
!!$          varname_spec = 'ss3_m'
!!$          cf_spec_stnd = 'seasalt_dry_aerosol_mode3_mass'
!!$          cf_spec_long = 'SS3_m'
!!$        case ( 'DUST2_N', 'dust2_n' )
!!$          varname_spec = 'dust2_n'
!!$          cf_spec_stnd = 'dust_dry_aerosol_mode2_number'
!!$          cf_spec_long = 'DUST2_n'
!!$        case ( 'DUST2_M', 'dust2_m' )
!!$          varname_spec = 'dust2_m'
!!$          cf_spec_stnd = 'dust_dry_aerosol_mode2_madust'
!!$          cf_spec_long = 'DUST2_m'
!!$        case ( 'DUST3_N', 'dust3_n' )
!!$          varname_spec = 'dust3_n'
!!$          cf_spec_stnd = 'dust_dry_aerosol_mode3_number'
!!$          cf_spec_long = 'DUST3_n'
!!$        case ( 'DUST3_M', 'dust3_m' )
!!$          varname_spec = 'dust3_m'
!!$          cf_spec_stnd = 'dust_dry_aerosol_mode3_madust'
!!$          cf_spec_long = 'DUST3_m'
        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


#ifdef with_m7
      end if ! RF%lpmx(k)
#endif       


      ! 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, sp_dat, temper_dat
#ifdef tropomi
    use MeteoData,   only : gph_dat
    use toolbox,     only : ltropo_ifs, lvlpress
#endif
#ifdef with_m7
    use calc_pm,     only : PMx_Integrate_3d
#endif
    !
    ! !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, i, j
    real, allocatable   ::  lev(:)
    integer             ::  l
    type(TDate)         ::  t, t0
    real                ::  time
    integer             ::  k, itr, dsec
    integer             ::  k_comp, itr_comp
    integer             ::  ims, ime, jms, jme, lms, lme
    integer             ::  gimr, gjmr, glmr

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

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

    ! for multiple of timestep only ...
    dsec = idate_f(4)*3600 + idate_f(5)*60 + idate_f(6)
    
    if ( modulo(dsec,RF%dsec) /= 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

#ifdef with_m7
    ! get helping pressure field in 3d
    allocate( pres3d(i1:i2,j1:j2,lmr) )
    ! fill mid level pressure
    call FPressure( levi, sp_dat(region)%data(i1:i2,j1:j2,1), pres3d, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! next time record:
    RF%trec = RF%trec + 1
    if(isRoot.and.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)

#ifdef tropomi
       ! As and Bs interfaces
       call MDF_Put_Var( RF%ncid, RF%varid_hyai, levi%a(0:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Var( RF%ncid, RF%varid_hybi, levi%b(0:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       ! As and Bs mid-level (full level)
       call MDF_Put_Var( RF%ncid, RF%varid_hyam, levi%fa(1:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       call MDF_Put_Var( RF%ncid, RF%varid_hybm, levi%fb(1:levi%nlev) , status)
       IF_NOTOK_MDF(fid=RF%ncid)
#else
       ! 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)
#endif

    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

    ! copy of temperature field
    RF%data3d_t(:,:,:,RF%trec) = temper_dat(region)%data(i1:i2,j1:j2,1:lmr)

    ! orography: copy of lowest interface gph field.  gph in the model is in "m", at interfaces, and gph(1)=oro
    ! only once ...    
    if ( RF%trec == 1 ) then
       RF%data2d_hgt(:,:) = gph_dat(region)%data(i1:i2,j1:j2,1)
    end if
    ! compute highest tropopause layer index
    do i = i1, i2
       do j = j1, j2
          RF%data2d_ltropo(i,j,RF%trec) = ltropo_ifs(region,i,j,temper_dat(region)%data(i,j,1:lmr),lmr)
       end do
    end do

#endif    

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

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


#ifdef with_m7

       ! ---------------------
       ! particulate matter
       ! ---------------------
       if( RF%lpmx(k) ) then 

          allocate( pmx( i1:i2, j1:j2, 1:lmr ) ) ; pmx = 0.0

          call PMx_Integrate_3d( region, RF%sizepmx(k), pmx, status )
          IF_NOTOK_RETURN(status=1)


          rf%data3d(:,:,:, rf%trec, k) = pmx
          
!           call MDF_Put_Var( RF%ncid, RF%varid_tr(k),                                           &
!                             reshape( pmx(i1:i2,j1:j2,lms:lme), (/imr,jmr,lmr,1/) ), status &
!                             start=(/i1,j1,1,RF%trec/), count=(/imr,jmr,lmr,1/) )
          deallocate( pmx )

       else 

#endif 
          
          ! ---------
          ! 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/) )
#ifdef with_m7
                
                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
#else
                write(gol,*)"Not using m7 - did not expected to be here."; call goErr
                write(gol,*)"             - make pres3d available"; call goErr
                status=1; TRACEBACK; return
#endif 
                
             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( 'numb' )
                ! write slab of concentrations
                !     number(trace)                     pressure            #/gridcell             Pa*K*mol
                ! C = ------------- * molmass_air * ----------------   =  ------------- * kg/mol *-----------
                !          m(air)                   temperature*Rgas       kg/gridcell                K*J
#ifdef with_m7                
                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
#else
                write(gol,*)"Not using m7 - did not expected to be here."; call goErr
                write(gol,*)"             - make pres3d available"; call goErr
                status=1; TRACEBACK; return
#endif

             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/) )
#ifdef with_m7                
                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
#else
                write(gol,*)"Not using m7 - did not expected to be here."; call goErr
                write(gol,*)"             - make pres3d available"; call goErr
                status=1; TRACEBACK; return
#endif
                
             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)

             ! ---------
             ! NOy and others (M7)
             ! ---------
#ifdef with_m7
          case( iNOy, iSO4, iBC, iPOM, iSS, iDU )
#else
          case( iNOy )
#endif 

            listtrac(:) = -999
            select case( itr )
            case( iNOy ); numtrac = nNOyt; listtrac(1:nNOyt) = iNOyt
#ifdef with_m7
            case( iSO4 ); numtrac = nSO4t; listtrac(1:nSO4t) = iSO4t
            case( iBC  ); numtrac = nBCt ; listtrac(1:nBCt ) = iBCt
            case( iPOM ); numtrac = nPOMt; listtrac(1:nPOMt) = iPOMt
            case( iSS  ); numtrac = nSSt ; listtrac(1:nSSt ) = iSSt
            case( iDU  ); numtrac = nDUt ; listtrac(1:nDUt ) = iDUt
#endif               
            end select

            ! mole fraction  = sum of mole fractions of components

            ! storage for sum of components (distributed over levels):
            allocate( Compo_k(i1:i2,j1:j2,lmr) )

            ! 3d fields with all levels or local levels only:
            allocate( field_k(i1:i2,j1:j2,lmr) )

            ! loop over transported components:
            Compo_k = 0.0
            do k_comp = 1, numtrac
               ! no more components?? 
               if( listtrac(k_comp) < 0 ) exit
               ! global tracer index:
               itr_comp = listtrac(k_comp)
               ! check ...
               if ( itr_comp > ntracet ) then
                  write (gol,'("index of NOy component does not represent a transported tracer : ",i3)') itr_comp; call goErr
                  TRACEBACK; status=1; return
               end if


               ! ----------------------------------------------------
               ! distinguish between mixing ratios and concentrations
               ! ----------------------------------------------------
               select case( RF%varid_type(k) )
               case( 'conc' )
                  ! calculate concentrations
                  !     m(trace)              pressure      xm(trace)
                  ! C = -------- * fscale * -----------  *  ---------
                  !      m(air)             temperature       Rgas
#ifdef with_m7                  
                  field_k = mass_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr_comp) / &
                            m_dat(region)%data(i1:i2,j1:j2,lms:lme)         * &
                            xmair * 1.E-03 * pres3d(i1:i2,j1:j2,1:lmr)      / &
                            temper_dat(region)%data(i1:i2,j1:j2,1:lmr)      / Rgas
#else
                  write(gol,*)"Not using m7 - did not expected to be here."; call goErr
                  write(gol,*)"             - make pres3d available"; call goErr
                  status=1; TRACEBACK; return
#endif 
               case( 'mixr' )
                  !     m(trace)
                  ! X = -------- * fscale 
                  !      m(air)
                  field_k = mass_dat(region)%rm(i1:i2,j1:j2,1:lmr,itr_comp) / &
                            m_dat(region)%data(i1:i2,j1:j2,lms:lme)         * &
                            fscale(itr_comp)
               case default 
                  write (gol,'("no such unit type",a)') RF%varid_type(k); call goErr
                  TRACEBACK; status=1; return
               end select

               ! add contribution of this component:
               Compo_k = Compo_k + field_k
            end do

            ! write slab of volume mixing ratio's:
!             call MDF_Put_Var( RF%ncid, RF%varid_tr(k), &
!                               reshape( Compo_k, (/imr,jmr,lmr,1/) ), &
!                               status, start=(/i1,j1,lms,RF%trec/), count=(/imr,jmr,lmr,1/) )
!             IF_NOTOK_MDF(fid=RF%ncid)

            rf%data3d(:,:,:, rf%trec, k) = Compo_k

            ! clear
            deallocate( Compo_k )
            deallocate( field_k )

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

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

         end select


#ifdef with_m7
      endif
#endif

   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)

      ! temperature (3d)
      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)      

#ifdef tropomi      

      if ( isRoot ) then
         write (gol,'(a,2i4)') 'PDUMP - writing fields T, hgt, ltropo, no2, so2, hcho;   trec, n_rec ', RF%trec, rf%n_rec
         call goPr
      end if

      ! temperature (3d)
      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)      

      ! surface altitude (2d)
      call MDF_Put_Var( RF%ncid, RF%varid_hgt, RF%data2d_hgt(:,:), status, start=(/i1,j1/), count=(/imr,jmr/) )
      IF_NOTOK_MDF(fid=RF%ncid)      

      ! highest tropopause level (2d)
      call MDF_Put_Var( RF%ncid, RF%varid_ltropo, RF%data2d_ltropo(:,:,:), status, start=(/i1,j1,1/), count=(/imr,jmr,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
   !----------------
#ifdef with_m7
   deallocate(pres3d)
#endif
   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 )
    deallocate(rf%data3d_t)

#ifdef tropomi
    deallocate(rf%data2d_hgt)
    deallocate(rf%data2d_ltropo)
#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
#ifdef with_m7
    RF%laod    = .false.
    RF%wavel   = -1.0
#endif
    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)

#ifdef with_m7

       
       ! ---------------------------
       ! check for AOD
       ! ---------------------------
       if( strlowercase(trname(1:3)) == 'aod' ) then 

          RF%laod(RF%ntr) = .true. 
          RF%itr (RF%ntr) = -1
          ! paste size to real 
          read(trname(5:len_trim(trname)), * )  RF%wavel(RF%ntr)

       else 

#endif

          ! convert to tm5 name:
          select case ( trim(strlowercase(trname)) )
          case ( 'hcho'        ) ; tmname = 'CH2O'
          case ( 'rn', 'radon' ) ; tmname = 'Rn222'
          case ( 'pb', 'lead'  ) ; tmname = 'Pb210'
          case default           ; tmname = trname
          end select

          ! NOy is a special ...
          select case ( trim(strlowercase(tmname)) )
          case ( 'NOy' )
             ! defined as ntrace+1
             RF%itr(RF%ntr) = iNOy
             write (gol,'("    *  ",a10)') trim(trname); call goPr
          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
#ifdef with_m7
       end if ! aod
#endif

       ! 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

#ifdef with_m7
       if( RF%laod(k) ) then 

          ! get diameter 
          write(cwavel,'(I3)') RF%wavel(k) 

          ! Aerosol Optical Depth (AOD): 
          varname_spec     = 'AOD@'//trim(cwavel)
          cf_spec_stnd     = 'AOD at '//trim(cwavel)//'nm'
          cf_spec_long     = 'aerosol optical depth at '//trim(cwavel)//' nanometer'
          cf_enti_stnd     = 'aerosol_optical_depth'
          cf_enti_unit     = '1'
          cf_enti_long     = 'aerosol optical depth'

       else 
#endif       


          ! 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 ( 'CO', 'co'   )
             varname_spec = 'co'
             cf_spec_stnd = 'carbon_monoxide'
             cf_spec_long = 'CO'
          case ( 'O3', 'o3'   )
             varname_spec = 'o3'
             cf_spec_stnd = 'ozone'
             cf_spec_long = 'O3'
          case ( 'O3s', 'o3s'   )
             varname_spec = 'o3s'
             cf_spec_stnd = 'ozone_from_stratosphere'
             cf_spec_long = 'O3s'
          case ( 'NO', 'no'  )
             varname_spec = 'no'
             cf_spec_stnd = 'nitrogen_monoxide'
             cf_spec_long = 'NO'
          case ( 'NO2', 'no2'  )
             varname_spec = 'no2'
             cf_spec_stnd = 'nitrogen_dioxide'
             cf_spec_long = 'NO2'
          case ( 'NOy', 'noy'  )
             varname_spec = 'noy'
             cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen'
             cf_spec_long = 'NOy'
             comment      = 'NOy = NOx + HNO3 + PAN + org.ntr., '// &
                  'with NOx = NO + NO2 + NO3 + HNO4 + N2O5'
          case ( 'CH2O', 'ch2o', 'CHOH', 'choh' )
             varname_spec = 'ch2o'
             cf_spec_stnd = 'formaldehyde'
             cf_spec_long = 'CH2O'
          case ( 'SO2', 'so2'  )
             varname_spec = 'so2'
             cf_spec_stnd = 'sulfur_dioxide'
             cf_spec_long = 'SO2'
          case ( 'CH4', 'ch4'  )
             varname_spec = 'ch4'
             cf_spec_stnd = 'methane'
             cf_spec_long = 'CH4'
          case ( 'OH', 'oh'  )
             varname_spec = 'oh'
             cf_spec_stnd = 'hydroxyl_radical'
             cf_spec_long = 'OH'
          case ( 'H2O2', 'h2o2'  )
             varname_spec = 'h2o2'
             cf_spec_stnd = 'hydrogen_peroxide'
             cf_spec_long = 'H2O2'
          case ( 'HNO3', 'hno3'  )
             varname_spec = 'hno3'
             cf_spec_stnd = 'nitric_acid'
             cf_spec_long = 'HNO3'
          case ( 'NH3', 'nh3'  )
             varname_spec = 'nh3'
             cf_spec_stnd = 'ammonia'
             cf_spec_long = 'NH3'
          case ( 'NH4', 'nh4'  )
             varname_spec = 'nh4'
             cf_spec_stnd = 'ammonium'
             cf_spec_long = 'NH4'
          case ( 'ORGNTR','orgntr'  )
             varname_spec = 'orgntr'
             cf_spec_stnd = 'organic_nitrate'
             cf_spec_long = 'ORGNTR'
          case ( 'PAN', 'pan'  )
             varname_spec = 'pan'
             cf_spec_stnd = 'peroxyacetyl_nitrate'
             cf_spec_long = 'PAN'
          case ( 'Rn', 'rn', 'Radon', 'radon'  )
             varname_spec = 'rn'
             cf_spec_stnd = 'radon'
             cf_spec_long = 'Rn'
          case ( 'Pb', 'pb', 'Lead', 'lead'  )
             varname_spec = 'pb'
             cf_spec_stnd = 'lead'
             cf_spec_long = 'Pb'
          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

#ifdef with_m7
       end if ! RF%laod(k)
#endif       

       ! 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
    integer(kind=8)   :: itau
    integer           :: 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)


!!$#ifdef with_m7
!!$
!!$                ! ---------------------
!!$                ! AOD
!!$                ! ---------------------
!!$                if( RF%laod(k) ) then 
!!$
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  allocate( ....... ) )
!!$
!!$                  call PMx_Integrate_3d( region, RF%sizepmx(k), pmx, status )
!!$                  IF_NOTOK_RETURN(status=1)
!!$                  
!!$                  ! root only:
!!$                  if ( myid == root ) then
!!$                    
!!$                    status = pnf90_put_var( RF%ncid, RF%varid_tr(k), &
!!$                         reshape( pmx(ims:ime,jms:jme,lms:lme), (/imr,jmr,lmr,1/) ), &
!!$                         start=(/1,1,1,RF%trec/), count=(/imr,jmr,lmr,1/) )
!!$                    
!!$                  end if
!!$                  
!!$                  deallocate( ............. )
!!$                  
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                  ! .............. PUT HERE THE CODE FOR AOD ..................! 
!!$                else 
!!$
!!$#endif 

                ! 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


#ifdef with_budgets
  
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE ##5 : 2D dry and wet deposition fields
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_DEPS_Init
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine RF_DEPS_Init( RF, fdir, model, expid, filetype, region, &
                           idate_f, dhour, 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
    !
    ! !OUTPUT PARAMETERS:
    !
    type(TPdumpFile_DEPS), intent(out) :: RF
    integer,               intent(out) :: status
    !
    ! !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
    !
    ! !REVISION HISTORY: 
    !    1 Oct 2010 - Achim Strunk - retor -> pdump
    !    7 Aug 2012 - Ph. Le Sager - netcdf4 thru mdf
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

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

    character(len=256)    ::  trnames
    character(len=8)      ::  trname, tmname
    integer               ::  k, itr

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

    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

    logical               ::  with_wdep
    integer               ::  imr, jmr, i1, i2, j1, j2

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

    call goLabel(rname)

    ! -- store arguments, init var
    RF%dhour        = dhour
    RF%tracer_names = tracer_names
    RF%ntr          = 0
    trnames         = tracer_names

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

    ! Switch to default .false., requires an extra call to PDUMP_Files_Write2 in OUTPUT_PDUMP_DONE
    n_deps_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, mess='DPS_Init' )
    !n_deps_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, .true., 'DPS_Init' )
    

    if ( n_deps_rec == 0 ) then     ! degenerated case
       deps_apply = .false.
       status=0
       return
    end if

    ! -- tracer index for requested tracers:
    if ( len_trim(trnames) == 0 ) then
       deps_apply = .false.
       write (gol,'("WARNING - NO tracers selected for depositions output!")') ; call goPr
       write (gol,'("        - deps_apply set to False.")'                   ) ; call goPr
       status=0
       return
    else
       write (gol,'("selected tracers for depositions output:")'); call goPr
    end if

    do
       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)

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

       ! convert to tm5 name:
       select case ( trname )
       case ( 'HCHO'        ) ; tmname = 'CH2O'
       case ( 'Rn', 'Radon' ) ; tmname = 'Rn222'
       case ( 'Pb', 'Lead'  ) ; tmname = 'Pb210'
       case default           ; tmname = trname
       end select

       ! wet deposition ?
       with_wdep = .false.
       select case ( trname )
       case ( 'HNO3' ) ; with_wdep = .true.
       case ( 'NOy'  ) ; with_wdep = .true.
       case ( 'NH3'  ) ; with_wdep = .true.
       case ( 'NH4'  ) ; with_wdep = .true.
       case ( 'SO4'  ) ; with_wdep = .true.
       end select
       RF%with_wdep(RF%ntr) = with_wdep

       ! NOy is a special ...
       select case ( tmname )
       case ( 'NOy' )
          ! defined as ntrace+1
          RF%itr(RF%ntr) = iNOy
          write (gol,'("  ",i3,"  ",a10,"  (",a10,") ",f12.4,"; wdep : ",l1)') &
               -1,trim(trname), '*', -1.0, with_wdep; call goPr
       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," ; wdep : ",l1)') &
                     itr, trim(trname), trim(names(itr)), ra(itr), with_wdep; 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
    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


    ! allocate storage:
    allocate( RF%ddep_budget(imr,jmr,RF%ntr) ) ; RF%ddep_budget = 0.0
    allocate( RF%wdep_budget(imr,jmr,RF%ntr) ) ; RF%wdep_budget = 0.0

    ! store current time (when budgets are reset):
    RF%t0_budget = NewDate(time6=idate_f)

    ! o open file

    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)

#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'              , 'dry and wet deposition' , 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, 'time'   , n_deps_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, '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, 'accum', 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  , 'long_name', 'length of accumulated time interval', status)
    IF_NOTOK_MDF(fid=RF%ncid)
    call MDF_Put_Att( RF%ncid, varid  , 'units'    , 'second'                             , status)
    IF_NOTOK_MDF(fid=RF%ncid)
    RF%varid_accum = varid

    allocate( RF%time(n_deps_rec)   )
    allocate( RF%date(6,n_deps_rec) )
    allocate( RF%dt(n_deps_rec)     )

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

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

       ! ~~ dry deposition

       ! CF standard name for concentration/mixing ratio/column:
       cf_enti_stnd = 'surface_dry_deposition_mole_flux'
       cf_enti_unit = 'mole m-2 s-1'
       cf_enti_long = 'dry deposition of '

       ! start of dataset name:
       varname_enti = 'dry'

       ! no comment yet
       comment = ''

       ! standard names from CF conventions:
       select case ( RF%name_tr(k) )
       case ( 'CO', 'co'   )
          varname_spec = 'co'
          cf_spec_stnd = 'carbon_monoxide'
          cf_spec_long = 'CO'
       case ( 'O3', 'o3'   )
          varname_spec = 'o3'
          cf_spec_stnd = 'ozone'
          cf_spec_long = 'O3'
       case ( 'O3s', 'o3s'   )
          varname_spec = 'o3s'
          cf_spec_stnd = 'ozone_from_stratosphere'
          cf_spec_long = 'O3s'
       case ( 'NO', 'no'  )
          varname_spec = 'no'
          cf_spec_stnd = 'nitrogen_monoxide'
          cf_spec_long = 'NO'
       case ( 'NO2', 'no2'  )
          varname_spec = 'no2'
          cf_spec_stnd = 'nitrogen_dioxide'
          cf_spec_long = 'NO2'
       case ( 'NOy', 'noy'  )
          varname_spec = 'noy'
          cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen'
          cf_spec_long = 'NOy'
          comment      = 'NOy = NOx + HNO3 + PAN + org.ntr., '// &
               'with NOx = NO + NO2 + NO3 + HNO4 + N2O5'
       case ( 'CH2O', 'ch2o', 'CHOH', 'choh' )
          varname_spec = 'ch2o'
          cf_spec_stnd = 'formaldehyde'
          cf_spec_long = 'CH2O'
       case ( 'SO2', 'so2'  )
          varname_spec = 'so2'
          cf_spec_stnd = 'sulfur_dioxide'
          cf_spec_long = 'SO2'
       case ( 'CH4', 'ch4'  )
          varname_spec = 'ch4'
          cf_spec_stnd = 'methane'
          cf_spec_long = 'CH4'
       case ( 'OH', 'oh'  )
          varname_spec = 'oh'
          cf_spec_stnd = 'hydroxyl_radical'
          cf_spec_long = 'OH'
       case ( 'H2O2', 'h2o2'  )
          varname_spec = 'h2o2'
          cf_spec_stnd = 'hydrogen_peroxide'
          cf_spec_long = 'H2O2'
       case ( 'HNO3', 'hno3'  )
          varname_spec = 'hno3'
          cf_spec_stnd = 'nitric_acid'
          cf_spec_long = 'HNO3'
       case ( 'NH3', 'nh3'  )
          varname_spec = 'nh3'
          cf_spec_stnd = 'ammonia'
          cf_spec_long = 'NH3'
       case ( 'ORGNTR','orgntr'  )
          varname_spec = 'orgntr'
          cf_spec_stnd = 'organic_nitrate'
          cf_spec_long = 'ORGNTR'
       case ( 'NH4', 'nh4'  )
          varname_spec = 'nh4'
          cf_spec_stnd = 'ammonium'
          cf_spec_long = 'NH4'
       case ( 'PAN', 'pan'  )
          varname_spec = 'pan'
          cf_spec_stnd = 'peroxyacetyl_nitrate'
          cf_spec_long = 'PAN'
       case ( 'Rn', 'rn', 'Radon', 'radon'  )
          varname_spec = 'rn'
          cf_spec_stnd = 'radon'
          cf_spec_long = 'Rn'
       case ( 'Pb', 'pb', 'Lead', 'lead'  )
          varname_spec = 'pb'
          cf_spec_stnd = 'lead'
          cf_spec_long = 'Pb'
       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

       write (varname,'(a,"_",a)') trim(varname_enti), trim(varname_spec)

       ! define variable:
       call MDF_Def_Var( RF%ncid, trim(varname), 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)    

       ! total names:
       cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)
       cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_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 ) 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_ddep(k) = varid

       
       ! ~~ wet deposition

       if ( RF%with_wdep(k) ) then

          ! CF standard name for concentration/mixing ratio/column:
          cf_enti_stnd = 'surface_wet_deposition_mole_flux'
          cf_enti_unit = 'mole m-2 s-1'
          cf_enti_long = 'wet deposition of '

          ! start of dataset name:
          varname_enti = 'wet'

          ! by default no comment:
          comment = ''

          ! standard names from CF conventions:
          select case ( RF%name_tr(k) )
          case ( 'NOy', 'noy'  )
             varname_spec = 'noy'
             cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen'
             cf_spec_long = 'NOy'
             comment      = 'NOy = NOx + HNO3 + PAN + org.ntr., '// &
                  'with NOx = NO + NO2 + NO3 + HNO4 + N2O5'
          case ( 'HNO3', 'hno3'  )
             varname_spec = 'hno3'
             cf_spec_stnd = 'nitric_acid'
             cf_spec_long = 'HNO3'
          case ( 'NH3', 'nh3'  )
             varname_spec = 'nh3'
             cf_spec_stnd = 'ammonia'
             cf_spec_long = 'NH3'
          case ( 'NH4', 'nh4'  )
             varname_spec = 'nh4'
             cf_spec_stnd = 'ammonium'
             cf_spec_long = 'NH4'
          case ( 'SO2', 'so2'  )
             varname_spec = 'so2'
             cf_spec_stnd = 'sulfur_dioxide'
             cf_spec_long = 'SO2'
          case default
             write (gol,'("unsupported tracer name for CF standard name : ",a)') RF%name_tr(k); call goPr
             TRACEBACK; status=1; return
          end select

          write (varname,'(a,"_",a)') trim(varname_enti), trim(varname_spec)

          ! define variable:
          call MDF_Def_Var( RF%ncid, trim(varname), 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)    

          ! total names:
          cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)
          cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_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 ) 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_wdep(k) = varid

       end if

    end do

    allocate( RF%data2d_dry(i1:i2, j1:j2, n_deps_rec, RF%ntr) )
    allocate( RF%data2d_wet(i1:i2, j1:j2, n_deps_rec, RF%ntr) )

!    RF%data2d_dry = 0. 
!    RF%data2d_wet = 0. 

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

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_DEPS_Write
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_DEPS_Write( RF, region, idate_f, status )
    !
    ! !USES:
    !
    use GO,             only : TDate, NewDate, Set, iTotal, rTotal, operator(-), wrtgol
    use Grid,           only : AreaOper
    use MeteoData,      only : global_lli, levi, lli
#ifndef without_chemistry
    use ebischeme,      only : buddrydep_dat => buddep_dat
#endif
#ifndef without_wet_deposition
    use wet_deposition, only : buddep_dat
#endif
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_DEPS), 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_DEPS_Write'

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

    integer             ::  imr, jmr, lmr
    type(TDate)         ::  t, t0
    real                ::  time
    real                ::  dt_sec
    integer             ::  k, itr, i1, i2, j1, j2

    real, allocatable   ::  budget(:,:)
    real, allocatable   ::  budget_loc(:,:)
    real, allocatable   ::  depflux(:,:)
    integer             ::  icomp

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

    ! temporary storage:
    allocate( budget_loc(imr,jmr) )
    allocate( depflux   (imr,jmr) )

    ! next time record:
    RF%trec = RF%trec + 1
    if(okdebug)then
       write(gol,*) "RF_DEPS_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' )

    ! length of time interval:
    dt_sec = iTotal( t - RF%t0_budget, 'sec' )

    ! zero time interval ? routine should not have been called ...
    if ( dt_sec == 0 ) then
      write (gol,'("routine called after zero lenght time interval:")'); call goErr
      call wrtgol( '  t0_budget : ', RF%t0_budget ); call goErr
      call wrtgol( '  t         : ', t            ); call goErr
      !status=1
      TRACEBACK
    end if

    ! reset timer:
    call Set( RF%t0_budget, time6=idate_f )

    !---------------
    ! Write GRID
    !---------------
    if ( RF%trec == 1 ) then

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

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

    end if

    !---------------    
    ! FILL DIAGNOSTIC ARRAYS
    !---------------
    
    !--------------- time 
    
    rf%time(rf%trec)   = time
    rf%date(:,rf%trec) = real(idate_f)
    rf%dt(rf%trec)     = dt_sec

    !--------------- dry deposition 
    
    do k = 1, RF%ntr

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

      ! extract current budget
#ifndef without_chemistry
      if ( itr == iNOy ) then
        ! add contributions of all NOy components:
        budget_loc = 0.0
        do icomp = 1, nNOyt
          budget_loc = budget_loc + buddrydep_dat(region)%dry(:,:,iNOyt(icomp))
        end do
      else
        ! extract budget for requested tracer:
        budget_loc = buddrydep_dat(region)%dry(:,:,itr)
      end if
#else
      budget_loc = 0.0
#endif
            
      ! deposition flux ~ (current budget - previous budget)/dt
      depflux = ( budget_loc - RF%ddep_budget(:,:,k) ) / dt_sec       ! mole/s
      call AreaOper( lli(region), depflux, '/', 'm2', status )    ! mole/m2/s
      IF_NOTOK_RETURN(status=1)

      ! save current budget & store record
      RF%ddep_budget(:,:,k) = budget_loc
      rf%data2d_dry(:,:,RF%trec,k)= depflux

    end do

    !--------------- wet deposition 

    do k = 1, RF%ntr

      ! skip ?
      if ( .not. RF%with_wdep(k) ) cycle

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

      ! extract current budget
#ifndef without_wet_deposition
      if ( itr == iNOy ) then
        ! add contributions of all NOy components:
        budget_loc = 0.0
        do icomp = 1, nNOyt
          ! add wet depositions for large scale and convective precip; total column:
          budget_loc = budget_loc + sum(buddep_dat(region)%lsp(:,:,:,iNOyt(icomp)),3) + &
                                    sum(buddep_dat(region)% cp(:,:,:,iNOyt(icomp)),3)
        end do
      else
        ! extract budget for requested tracer;
        ! add wet depositions for large scale and convective precip; total column:
        budget_loc = sum(buddep_dat(region)%lsp(:,:,:,itr),3) + &
                     sum(buddep_dat(region)% cp(:,:,:,itr),3)
      end if
#else
      budget_loc = 0.0
#endif

      ! deposition flux ~ (current budget - previous budget)/dt
      depflux = ( budget_loc - RF%wdep_budget(:,:,k) ) / dt_sec    ! mole/s
      call AreaOper( lli(region), depflux, '/', 'm2', status )     ! mole/m2/s
      IF_NOTOK_RETURN(status=1)

      ! save current budget & store record
      RF%wdep_budget(:,:,k)       = budget_loc
      RF%data2d_wet(:,:,RF%trec,k)= depflux


    end do

    !----------------
    !    WRITE
    !----------------
    if ( RF%trec == n_deps_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)

       ! accumulation interval
       call MDF_Put_Var( RF%ncid, RF%varid_accum, rf%dt , status)
       IF_NOTOK_MDF(fid=RF%ncid)

       ! deposition flux
       do k = 1, RF%ntr

          call MDF_Put_Var( RF%ncid, RF%varid_ddep(k), rf%data2d_dry(:,:,:,k), status, start=(/i1,j1,1/) )
          IF_NOTOK_MDF(fid=RF%ncid)

          if ( .not. RF%with_wdep(k) ) cycle

          call MDF_Put_Var( RF%ncid, RF%varid_wdep(k), rf%data2d_wet(:,:,:,k), status, start=(/i1,j1,1/) )
          IF_NOTOK_MDF(fid=RF%ncid)

       end do

    end if

    
    !----------------
    ! DONE
    !----------------
    deallocate( budget_loc )
    deallocate( depflux    )
    call goLabel()
    status = 0

  END SUBROUTINE RF_DEPS_Write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RF_DEPS_Done
  !
  ! !DESCRIPTION: close file #5
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_DEPS_Done( RF, status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_DEPS), 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_DEPS_Done'

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

    call goLabel(rname)

    ! close file
    call MDF_Close( RF%ncid , status)
    IF_NOTOK_RETURN(status=1)
    
    ! clear
    deallocate( RF%ddep_budget )
    deallocate( RF%wdep_budget )
    deallocate( rf%time, rf%date, rf%dt, rf%data2d_dry, rf%data2d_wet )

    call goLabel() ; status = 0

  END SUBROUTINE RF_DEPS_Done
  !EOC



  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE #6 : deposition velocities
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  RF_DEPV_Init
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine RF_DEPV_Init( RF, fdir, model, expid, filetype, region, &
                           idate_f, dhour, 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
    !
    ! !OUTPUT PARAMETERS:
    !
    type(TPdumpFile_DEPV), 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)                    ::  dhour
    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_DEPV_Init'

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

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

    character(len=256)    ::  trnames
    character(len=8)      ::  trname, tmname
    integer               ::  k, itr

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

    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%dhour        = dhour
    RF%tracer_names = tracer_names
    RF%ntr          = 0
    trnames         = tracer_names

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

    n_depv_rec = GET_N_TIME_RECORDS( idate_f, dhour*3600, mess='DEPV_Init' )

    if ( n_depv_rec == 0 ) then  ! degenerated cases
       depv_apply = .false.
       status=0
       return
    end if

    ! tracer index for requested tracers
    if ( len_trim(trnames) == 0 ) then
       depv_apply = .false.
       write (gol,'("WARNING - NO tracers selected for depositions velocity output!")') ; call goPr
       write (gol,'("        - depv_apply set to False.")'                            ) ; call goPr
       status=0
       return
    else
       write (gol,'("selected tracers for deposition velocity output:")'); call goPr       
    end if

    
    do
       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)

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

       ! convert to tm5 name:
       select case ( trname )
       case ( 'HCHO'        ) ; tmname = 'CH2O'
       case ( 'Rn', 'Radon' ) ; tmname = 'Rn222'
       case ( 'Pb', 'Lead'  ) ; tmname = 'Pb210'
       case default           ; tmname = trname
       end select

       ! 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

       ! 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
    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'              , 'volume mixing ratios'   , 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, 'time'   , n_depv_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, '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_depv_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_depv_rec) )

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

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

       ! CF standard name for concentration/mixing ratio/column:
       cf_enti_stnd = 'surface_dry_deposition_velocity_due_to_turbulence'
       cf_enti_unit = 'mole m-2 s-1'
       cf_enti_long = 'dry deposition of '

       ! start of dataset name:
       varname_enti = 'ddepvel'

       ! by default no comment:
       comment = ''

       ! standard names from CF conventions:
       select case ( RF%name_tr(k) )
       case ( 'CO', 'co'   )
          varname_spec = 'co'
          cf_spec_stnd = 'carbon_monoxide'
          cf_spec_long = 'CO'
       case ( 'O3', 'o3'   )
          varname_spec = 'o3'
          cf_spec_stnd = 'ozone'
          cf_spec_long = 'O3'
       case ( 'O3s', 'o3s'   )
          varname_spec = 'o3s'
          cf_spec_stnd = 'ozone_from_stratosphere'
          cf_spec_long = 'O3s'
       case ( 'NO', 'no'  )
          varname_spec = 'no'
          cf_spec_stnd = 'nitrogen_monoxide'
          cf_spec_long = 'NO'
       case ( 'NO2', 'no2'  )
          varname_spec = 'no2'
          cf_spec_stnd = 'nitrogen_dioxide'
          cf_spec_long = 'NO2'
       case ( 'NOy', 'noy'  )
          varname_spec = 'noy'
          cf_spec_stnd = 'all_nitrogen_oxides_as_nitrogen'
          cf_spec_long = 'NOy'
          comment      = 'NOy = NOx + HNO3 + PAN + org.ntr., '// &
               'with NOx = NO + NO2 + NO3 + HNO4 + N2O5'
       case ( 'CH2O', 'ch2o', 'CHOH', 'choh' )
          varname_spec = 'ch2o'
          cf_spec_stnd = 'formaldehyde'
          cf_spec_long = 'CH2O'
       case ( 'SO2', 'so2'  )
          varname_spec = 'so2'
          cf_spec_stnd = 'sulfur_dioxide'
          cf_spec_long = 'SO2'
       case ( 'CH4', 'ch4'  )
          varname_spec = 'ch4'
          cf_spec_stnd = 'methane'
          cf_spec_long = 'CH4'
       case ( 'OH', 'oh'  )
          varname_spec = 'oh'
          cf_spec_stnd = 'hydroxyl_radical'
          cf_spec_long = 'OH'
       case ( 'H2O2', 'h2o2'  )
          varname_spec = 'h2o2'
          cf_spec_stnd = 'hydrogen_peroxide'
          cf_spec_long = 'H2O2'
       case ( 'HNO3', 'hno3'  )
          varname_spec = 'hno3'
          cf_spec_stnd = 'nitric_acid'
          cf_spec_long = 'HNO3'
       case ( 'PAN', 'pan'  )
          varname_spec = 'pan'
          cf_spec_stnd = 'peroxyacetyl_nitrate'
          cf_spec_long = 'PAN'
       case ( 'Rn', 'rn', 'Radon', 'radon'  )
          varname_spec = 'rn'
          cf_spec_stnd = 'radon'
          cf_spec_long = 'Rn'
       case ( 'Pb', 'pb', 'Lead', 'lead'  )
          varname_spec = 'pb'
          cf_spec_stnd = 'lead'
          cf_spec_long = 'Pb'
       case ( 'NH3', 'nh3'  )
          varname_spec = 'nh3'
          cf_spec_stnd = 'ammonia'
          cf_spec_long = 'NH3'
       case ( 'NH4', 'nh4'  )
          varname_spec = 'nh4'
          cf_spec_stnd = 'ammonium'
          cf_spec_long = 'NH4'
       case default
          write (gol,'("unsupported tracer name for CF standard name : ",a)') RF%name_tr(k); call goPr
          TRACEBACK; status=1; return
       end select

       write (varname,'(a,"_",a)') trim(varname_enti), trim(varname_spec)    

       write (gol,'("    varname : ",a)') trim(varname); call goPr

       ! define variable:
       call MDF_Def_Var( RF%ncid, trim(varname), 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)    

       ! total names:
       cf_name_stnd = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)
       cf_name_long = trim(cf_enti_long)//' of '//trim(cf_spec_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)
       call MDF_Put_Att( RF%ncid   , varid, 'moleweight_tracer', ra(itr)*1e3       , status)
       IF_NOTOK_MDF(fid=RF%ncid)
       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

    allocate( rf%data2d(i1:i2, j1:j2, n_depv_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_DEPV_Init
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  RF_DEPV_Write
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_DEPV_Write( RF, region, idate_f, status )
    !
    ! !USES:
    !
    use GO,             only : TDate, NewDate, Set, iTotal, rTotal, operator(-), wrtgol
    use Grid,           only : AreaOper
    use MeteoData,      only : global_lli
#ifndef without_dry_deposition
    use dry_deposition, only : vd
#endif
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_DEPV), 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_DEPV_Write'

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

    integer             ::  imr, jmr
    type(TDate)         ::  t, t0
    real                ::  time
    integer             ::  k, itr, i1, i2, j1, j2
    real, allocatable   ::  depvel(:,:)

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

    ! next time record:
    RF%trec = RF%trec + 1
    if(okdebug)then
       write(gol,*) "RF_DEPV_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 : Dimensions    
    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)

    end if

    !-------- FILL DIAGNOSTIC ARRAYS

    rf%time(rf%trec) = time
    rf%date(:,rf%trec) = real(idate_f)
    
    ! loop over tracers to be written:
    do k = 1, RF%ntr

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

#ifndef without_dry_deposition
       rf%data2d(:,:,RF%trec,k) = vd(region,itr)%surf   ! deposition velocity
#else
       rf%data2d(:,:,RF%trec,k) = 0.0
#endif
    end do

    
    !-------- WRITE
    if ( RF%trec == n_depv_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)

       ! loop over tracers to be written:
       do k = 1, RF%ntr
          call MDF_Put_Var( RF%ncid, RF%varid_tr(k), rf%data2d(:,:,:,k), status, start=(/i1,j1,1/))
          IF_NOTOK_MDF(fid=RF%ncid)
       end do
       
    end if
    

    call goLabel()
    status = 0

  END SUBROUTINE RF_DEPV_Write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  RF_DEPV_Done
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RF_DEPV_Done( RF, status )
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    type(TPdumpFile_DEPV), 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_DEPV_Done'

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

    call goLabel(rname)

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

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

    call goLabel() ; status = 0

  END SUBROUTINE RF_DEPV_Done
  !EOC
  
#endif
  
  !--------------------------------------------------------------------------
  !                    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