ecearth3/sources/tm5mp/proj/output/user_output_cf.F90

!#################################################################
!
! CF-conform HDF Output
!
! NetCDF files following CF conventions.
!
! Based on CF conform output see:
!     website
!     http://cf-pcmdi.llnl.gov/
!
!    NetCDF Climate and Forecast (CF) Metadata Convention
!   
!
!   Files should now follow the CF conventions, which:
!
!    - longitudes from [0,360] ?
!    - levels from surface to top
!    - time from 2001-01-01 00:00
!
! SAMPLE RCFILE
!
!    output.cf                        :  T
!    output.cf.dataset.author         :  A.  Slave (KNMI)
!    output.cf.dataset.institution    :  KNMI, De Bilt, The Netherlands
!    output.cf.dataset.version        :  GEMS run
!    output.cf.fname.model            :  TM5
!    output.cf.fname.expid            :  V1
!
!    output.cf.griddef.apply          :  F
!
!    output.cf.tp.apply               :  T
!
!    output.cf.vmr.n                  :  3
!    output.cf.vmr.001.apply          :  T
!    output.cf.vmr.001.fname          :  vmr3
!    output.cf.vmr.001.dhour          :  3
!    output.cf.vmr.001.tracers        :  SO2 NOy CH4 OH HNO3 PAN H2O2 Radon Lead
!    output.cf.vmr.002.apply          :  T
!
!    output.cf.vmr.002.fname          :  vmr1
!    output.cf.vmr.002.dhour          :  1
!    output.cf.vmr.002.tracers        :  O3 O3s CO NO2 NO CH2O
!    output.cf.vmr.003.apply          :  F
!
!    output.cf.vmr.003.fname          :  vmra
!    output.cf.vmr.003.dhour          :  3
!    output.cf.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.cf.depositions.apply      :  F
!    output.cf.depositions.dhour      :  3
!    output.cf.depositions.tracers    :  O3 HNO3 NO NO2 H2O2 CH2O PAN CO NH3 NH4 SO2 NOy
!
!    output.cf.depvels.apply          :  F
!    output.cf.depvels.dhour          :  3
!    output.cf.depvels.tracers        :  O3 HNO3 NO NO2 H2O2 CH2O PAN CO NH3 NH4 SO2
!    
!
!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################

module user_output_cf

  use GO          , only : gol, goPr, goErr, goLabel
  use GO          , only : TDate
  use dims        , only : nregions
  use chem_param  , only : ntrace
  use chem_param  , only : iNOx, iHNO3, iNO3_a, iPAN, iOrgNtr
#ifdef with_tendencies
  use tm5_tendency, only : plc_ntr, plc_npr
#endif   

  use file_hdf,    only : THdfFile, TSds,setdim
  
  implicit none


  ! --- in/out --------------------------
  
  private
  
  public  ::  Output_CF_Init, Output_CF_Done, Output_CF_Step
  
  
  ! --- const ---------------------------
  
  character(len=*), parameter  ::  mname = 'user_output.cf'
  
  ! reference time:
  integer, parameter            ::  time_reftime6(6) = (/2000,01,01,00,00,00/)
  character(len=*), parameter   ::  time_units = 'days since 2000-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 
  !
  integer, parameter  ::  iNOy = ntrace + 1    
  integer, parameter  ::  nNOyt = 5
  integer, parameter  ::  iNOyt(nNOyt) = (/ iNOx, iHNO3, iNO3_a, iPAN, iOrgNtr /)
  
  ! --- types ---------------------------
  
  type TCF_File_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 TCF_File_GridDef
  
  type TCF_File_TP
    integer             ::  trec
    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
  end type TCF_File_TP
  
  type TCF_File_VMR
    integer             ::  trec
    integer             ::  dhour
    character(len=256)  ::  tracer_names
    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             ::  ntr
    integer             ::  itr(ntrace)
    character(len=8)    ::  name_tr(ntrace)
    integer             ::  varid_tr(ntrace)
  end type TCF_File_VMR
   
  type TCF_File_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
  end type TCF_File_DEPS
  
  type TCF_File_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)
  end type TCF_File_DEPV

#ifdef with_tendencies
  type TCF_File_TEND
    integer             ::  trec
    integer             ::  dhour
    character(len=256)  ::  tracer_names
    character(len=256)  ::  proces_names
    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             ::  ntr
    integer             ::  itr(plc_ntr)
    character(len=8)    ::  name_tr(plc_ntr)
    integer             ::  npr
    integer             ::  ipr(plc_npr)
    character(len=8)    ::  name_pr(plc_npr)
    integer             ::  varid_tend(plc_ntr,plc_npr)
  end type TCF_File_TEND
#endif   
  type(THdfFile)        :: GridDef_hdf(nregions)

 type TCF_output_var
    integer             ::  ivar
    character(len=128)  ::  var_name
    character(len=32)   ::  param_name
    character(len=32)   ::  param_name_default
    character(len=256)  ::  stand_name
    character(len=256)  ::  stand_name_default
    character(len=512)  ::  long_name
    character(len=512)  ::  long_name_default
    character(len=256)  ::  units
    character(len=256)  ::  units_default
  end type TCF_output_var

  ! --- var -----------------------------
  
  integer                         ::  curr_day(nregions,3)
  
  character(len=32)               ::  fname_model
  character(len=6)                ::  fname_expid
  character(len=32)               ::  fname_grid(nregions)
  character(len=256)              ::  dataset_author, dataset_institution, dataset_version
  
  logical                                      ::  griddef_apply
  type(TCF_File_GridDef), save                 ::  CFGridDef(nregions)

  integer                                      ::  nvar_dims

  integer,parameter                            ::  CF_grid_nvar=10
  type(TCF_output_var), save                   ::  CF_grid_var(CF_grid_nvar)
  type(TSds), save                             ::  sds_grid(CF_grid_nvar) 

  
  logical, allocatable                         ::  vmr_apply(:)
  type(TCF_output_var), save                   ::  CF_vmr_var(ntrace)
  type(TSds), allocatable, save                ::  sds_vmr_inst(:,:,:) 
  
  logical                                      ::  tp_apply
  ! currently only surface pressure; to be extended
  type(TCF_output_var), save                   ::  CF_TP_var(1) 

  integer, save                                ::  nvmr
  logical, allocatable                         ::  vmr_apply_average(:)
  character(len=16), allocatable               ::  vmr_fname(:)
  integer, allocatable                         ::  vmr_nvars(:)
  integer, allocatable                         ::  vmr_dhour(:)
  character(len=256), allocatable              ::  vmr_tracer_names(:)
  type(TCF_File_VMR), allocatable, save        ::  CFVMR(:,:)
  type(THdfFile), allocatable, save            ::  CFVMR_hdf_inst(:,:)
  type(THdfFile), allocatable, save            ::  CFVMR_hdf_average(:,:)
  type(TCF_File_TP), save                      ::  CFTP(nregions)
  type(TSds), allocatable, save                ::  sds_vmr_average(:,:,:) 
   
  logical                                      ::  deps_apply
  character(len=16)                            ::  deps_fname
  integer                                      ::  deps_dhour
  character(len=256)                           ::  deps_tracer_names
  type(TCF_File_DEPS), save                    ::  CFDEPS(nregions)
  
  logical                                      ::  depv_apply
  character(len=16)                            ::  depv_fname
  integer                                      ::  depv_dhour
  character(len=256)                           ::  depv_tracer_names
  type(TCF_File_DEPV), save                    ::  CFDEPV(nregions)


contains


  ! ********************************************************************
  ! ***
  ! *** init/step/done cf output
  ! *** and the naming
  ! ***
  ! ********************************************************************
  
  ! ***
  
  ! ********************************************************************
  ! ***
  ! *** subroutine Output_CF_naming reads in the rc-file that 
  ! *** provides the user-defined naming conventions for the variables that 
  ! *** are requested to be written out.
  ! *** The attributes currently contain the parameter name, default name, long name, 
  ! *** and units.
  ! *** 
  ! *** All default settings can be overwritten by the data from the rc-file.
  ! ********************************************************************

  subroutine Output_CF_naming(rcF,status)
      use GO   , only : TrcFile, ReadRc
      use chem_param, only : ntrace, names, ra
   ! --- in/out ---------------------------------
    
    type(TrcFile), intent(inout)     ::  rcF
    integer, intent(out)             ::  status

   ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Output_CF_naming'
    
    ! --- local ------------------------------
    
    character(len=128)        ::  test_param_name
    integer                   ::  ivar
    
    character(len=32)     ::  varname, varname_conc, 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 -------------------------------
   !initialize standard names of keys...

    ivar=1
    CF_grid_var(ivar)%var_name='lon'
    CF_grid_var(ivar)%param_name_default='lon'
    CF_grid_var(ivar)%stand_name_default='longitude'
    CF_grid_var(ivar)%long_name_default='longitude'
    CF_grid_var(ivar)%units_default='degrees_east'
         
    ivar=2
    CF_grid_var(ivar)%var_name='lat'
    CF_grid_var(ivar)%param_name_default='lat'
    CF_grid_var(ivar)%stand_name_default='latitude'
    CF_grid_var(ivar)%long_name_default='latitude'
    CF_grid_var(ivar)%units_default='degrees_north'

    ivar=3
    CF_grid_var(ivar)%var_name='lev'
    CF_grid_var(ivar)%param_name_default='lev'
    CF_grid_var(ivar)%stand_name_default='level'
    CF_grid_var(ivar)%long_name_default='level'
    CF_grid_var(ivar)%units_default='level'

    ivar=4
    CF_grid_var(ivar)%var_name='time'
    CF_grid_var(ivar)%param_name_default='time'
    CF_grid_var(ivar)%stand_name_default='time'
    CF_grid_var(ivar)%long_name_default='time'
    CF_grid_var(ivar)%units_default='time'

    ivar=5
    CF_grid_var(ivar)%var_name='date'
    CF_grid_var(ivar)%param_name_default='date'
    CF_grid_var(ivar)%stand_name_default='date'
    CF_grid_var(ivar)%long_name_default='date'
    CF_grid_var(ivar)%units_default='date'

    ivar=6
    CF_grid_var(ivar)%var_name='area'
    CF_grid_var(ivar)%param_name_default='area'
    CF_grid_var(ivar)%stand_name_default='grid_cell_area'
    CF_grid_var(ivar)%long_name_default='grid-cell area'
    CF_grid_var(ivar)%units_default='m2'         

    ivar=7
    CF_grid_var(ivar)%var_name='a'
    CF_grid_var(ivar)%param_name_default='a'
    CF_grid_var(ivar)%stand_name_default='atmosphere_hybrid_sigma_coordinate'
    CF_grid_var(ivar)%long_name_default='hybrid sigma coordinate a coefficient'
    CF_grid_var(ivar)%units_default='1'
    
    ivar=8
    CF_grid_var(ivar)%var_name='b'
    CF_grid_var(ivar)%param_name_default='b'
    CF_grid_var(ivar)%stand_name_default='atmosphere_hybrid_sigma_coordinate'
    CF_grid_var(ivar)%long_name_default='hybrid sigma coordinate b coefficient'
    CF_grid_var(ivar)%units_default='1'
    
    ivar=9
    CF_grid_var(ivar)%var_name='a_bnds'
    CF_grid_var(ivar)%param_name_default='a_bnds'
    CF_grid_var(ivar)%stand_name_default='atmosphere_hybrid_sigma_coordinate'
    CF_grid_var(ivar)%long_name_default='hybrid sigma coordinate a coefficient for layer bounds'
    CF_grid_var(ivar)%units_default='1'
    
    ivar=10
    CF_grid_var(ivar)%var_name='b_bnds'
    CF_grid_var(ivar)%param_name_default='b_bnds'
    CF_grid_var(ivar)%stand_name_default='atmosphere_hybrid_sigma_coordinate'
    CF_grid_var(ivar)%long_name_default='hybrid sigma coordinate b coefficient for layer bounds'
    CF_grid_var(ivar)%units_default='1'
        
    !read names of keys from rc-file
    do ivar =1,CF_grid_nvar
      call ReadRc( rcF, 'output.cf.griddef.'//trim( CF_grid_var(ivar)%var_name)//'.param.name', CF_grid_var(ivar)%param_name, status,default=trim(CF_grid_var(ivar)%param_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.griddef.'//trim( CF_grid_var(ivar)%var_name)//'.stand.name', CF_grid_var(ivar)%stand_name, status,default=trim(CF_grid_var(ivar)%stand_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.griddef.'//trim( CF_grid_var(ivar)%var_name)//'.long.name', CF_grid_var(ivar)%long_name  , status,default=trim(CF_grid_var(ivar)%long_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.griddef.'//trim( CF_grid_var(ivar)%var_name)//'.units'    , CF_grid_var(ivar)%units      , status,default=trim(CF_grid_var(ivar)%units_default) )
      IF_ERROR_RETURN(status=1)
    enddo
    

      ! loop over all available tracers to define
      ! CF standard name for species:

      cf_enti_stnd = 'mole_fraction' ; cf_enti_unit = 'mole mole-1' ; cf_enti_long = 'volume mixing ratio'
      
      ! CF standard name for medium:
      cf_medium_stnd = 'in_air' ; cf_medium_long = 'in humid air'

      ! begin of variable name:
      varname_conc = 'vmr'
      
      ! no comment yet
      comment = ''
      
      ! standard names from CF conventions:

    do ivar = 1,ntrace
      select case ( trim(names(ivar)) )
        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 = 'nitrogen_oxides'
          cf_spec_long = 'NOy'
          comment      = 'NOy = NOx + HNO3 + PAN + org.ntr., '// &
                           'with NOx = NO + NO2 + NO3 + HNO4 + N2O5'
        case ( 'NOx', 'nox'  )
          varname_spec = 'nox'
          cf_spec_stnd = 'nitrogen_oxides'
          cf_spec_long = 'NOx'
          comment      = '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 ( 'SO4', 'so4'  )
          varname_spec = 'so4'
          cf_spec_stnd = 'sulfate_as_sulfate_dry_aerosol'
          cf_spec_long = 'SO4'
        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 ( 'NO3_A', 'no3_a'  )
          varname_spec = 'no3'
          cf_spec_stnd = 'nitrate_as_nitrate_dry_aerosol'
          cf_spec_long = 'NO3'
        case ( 'BC', '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
          varname_spec = trim(names(ivar))
          cf_spec_stnd = trim(names(ivar))
          cf_spec_long = trim(names(ivar))

          write (gol,'("Use dummy naming for tracer: ",a)') names(ivar); call goPr
          write (gol,'("Use standard name")'); call goPr
!          TRACEBACK; status=1; return
      end select
        
      write (varname,'(a,"_",a)') trim(varname_conc), trim(varname_spec)
      
      
      ! construct total names:
      CF_vmr_var(ivar)%var_name           = trim(names(ivar))  ! e.g. O3, NO2
      CF_vmr_var(ivar)%param_name_default = trim(names(ivar))  ! the name of the output-variable. Currently the same, but may be different, e.g. vmr_o3, vmr_no2
      CF_vmr_var(ivar)%stand_name_default = trim(cf_enti_stnd)//'_of_'//trim(cf_spec_stnd)//'_'//trim(cf_medium_stnd)
      CF_vmr_var(ivar)%long_name_default  = trim(cf_enti_long)//' of '//trim(cf_spec_long)//' '//trim(cf_medium_long)
      CF_vmr_var(ivar)%units_default       = trim(cf_enti_unit)
    
      !overwrite these names of keys from rc-file, if necessary
      call ReadRc( rcF, 'output.cf.vmr.'//trim( CF_vmr_var(ivar)%var_name)//'.param.name', CF_vmr_var(ivar)%param_name, status,default=trim(CF_vmr_var(ivar)%param_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.vmr.'//trim( CF_vmr_var(ivar)%var_name)//'.stand.name', CF_vmr_var(ivar)%stand_name, status,default=trim(CF_vmr_var(ivar)%stand_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.vmr.'//trim( CF_vmr_var(ivar)%var_name)//'.long.name', CF_vmr_var(ivar)%long_name  , status,default=trim(CF_vmr_var(ivar)%long_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.vmr.'//trim( CF_vmr_var(ivar)%var_name)//'.units'    , CF_vmr_var(ivar)%units      , status,default=trim(CF_vmr_var(ivar)%units_default) )
      IF_ERROR_RETURN(status=1)
    enddo
      
      !
      !Define naming of selected meteo-variables:
      !
      ivar = 1   ! surface pressure
      
      CF_TP_var(ivar)%var_name           = 'ps'
      CF_TP_var(ivar)%param_name_default = 'ps'
      CF_TP_var(ivar)%stand_name_default = 'surface_air_pressure'
      CF_TP_var(ivar)%long_name_default  = 'surface pressure'
      CF_TP_var(ivar)%units_default       = 'Pa'
    
      !overwrite these names of keys from rc-file, if necessary
      call ReadRc( rcF, 'output.cf.TP.'//trim( CF_TP_var(ivar)%var_name)//'.param.name', CF_TP_var(ivar)%param_name, status,default=trim(CF_TP_var(ivar)%param_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.TP.'//trim( CF_TP_var(ivar)%var_name)//'.stand.name', CF_TP_var(ivar)%stand_name, status,default=trim(CF_TP_var(ivar)%stand_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.TP.'//trim( CF_TP_var(ivar)%var_name)//'.long.name', CF_TP_var(ivar)%long_name  , status,default=trim(CF_TP_var(ivar)%long_name_default) )
      IF_ERROR_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.TP.'//trim( CF_TP_var(ivar)%var_name)//'.units'    , CF_TP_var(ivar)%units      , status,default=trim(CF_TP_var(ivar)%units_default) )
      IF_ERROR_RETURN(status=1)
    
    
  end subroutine Output_CF_naming

  ! ***
  
  subroutine Output_CF_Init( rcF, dhour_min, status )
  
    use GO   , only : TrcFile, ReadRc
    use meteo, only : lli
    use partools   , only : Par_Barrier
    
    ! --- in/out ---------------------------------
    
    type(TrcFile), intent(inout)     ::  rcF
    integer, intent(out)             ::  dhour_min
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Output_CF_Init'
    
    ! --- local ------------------------------
    
    integer                   ::  region
    integer                   ::  itf
    character(len=64)         ::  key
    character(len=3)          ::  nr
    integer                   ::  ivmr,ivar,vmr_nvars_max
    character(len=256)        ::  test_param_name

    ! --- begin -------------------------------
    
    call goLabel(rname)

    ! Read in the dataset keys:
    call ReadRc( rcF, 'output.cf.dataset.author'     , dataset_author     , status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.cf.dataset.institution', dataset_institution, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.cf.dataset.version'    , dataset_version    , status )
    IF_NOTOK_RETURN(status=1)
    
    ! Read in the filename keys:
    call ReadRc( rcF, 'output.cf.fname.model', fname_model, status )
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'output.cf.fname.expid', fname_expid, status )
    IF_NOTOK_RETURN(status=1)
    
    ! Read in the 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.cf.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
    
    !
    ! Initialize CF-compliant naming conventions: Read in the naming from an rc-file 
    !
    call Output_CF_naming(rcF,status)
    IF_ERROR_RETURN(status=1)
    !
    !
    !

    
    ! Write-out a griddef file ?
    call ReadRc( rcF, 'output.cf.griddef.apply', griddef_apply, status )
    IF_NOTOK_RETURN(status=1)
    
    ! Write out a temperature, pressure, etc file ? (not yet implemented!)
    call ReadRc( rcF, 'output.cf.tp.apply', tp_apply, status )
    IF_NOTOK_RETURN(status=1)
    
    ! VMR files:
    ! read number of vmr files to be written:
    call ReadRc( rcF, 'output.cf.vmr.n', nvmr, status )
    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

    ! write any vmr files (data-files containing instantaneous info on species vol.mix.ratios)
    if ( nvmr > 0 ) then
      ! storage:
      allocate( vmr_apply(nvmr)         )
      allocate( vmr_apply_average(nvmr) )
      allocate( vmr_fname(nvmr)         )
      allocate( vmr_nvars(nvmr)         )
      allocate( vmr_dhour(nvmr)         )
      allocate( vmr_tracer_names(nvmr)  )
      allocate( CFVMR(nregions,nvmr)    )
      ! loop over vmr files:
      do ivmr = 1, nvmr
        ! number in rc keys:
        write (nr,'(i3.3)') ivmr
        ! apply ?
        call ReadRc( rcF, 'output.cf.vmr.'//nr//'.apply', vmr_apply(ivmr), status )
        IF_NOTOK_RETURN(status=1)
        ! proceed ?
        if ( vmr_apply(ivmr) ) then
          ! first part of filename:
          call ReadRc( rcF, 'output.cf.vmr.'//nr//'.fname', vmr_fname(ivmr), status )
          IF_NOTOK_RETURN(status=1)
          ! How many parameters are defined that need to be written out?
          call ReadRc( rcF, 'output.cf.vmr.'//nr//'.nvars', vmr_nvars(ivmr), status )
          IF_NOTOK_RETURN(status=1)

          ! time resolution:
          call ReadRc( rcF, 'output.cf.vmr.'//nr//'.dhour', vmr_dhour(ivmr), status )
          IF_NOTOK_RETURN(status=1)

          ! the names of tracers to be written:
          call ReadRc( rcF, 'output.cf.vmr.'//nr//'.tracers', vmr_tracer_names(ivmr), status )
          IF_NOTOK_RETURN(status=1)
        end if  ! apply ?
      end do  ! vmr numbers
      
      ! Find maximum number of variables to be saved:
      vmr_nvars_max = 0
      do ivmr=1,nvmr 
        vmr_nvars_max=max(vmr_nvars_max,vmr_nvars(ivmr))
      enddo

      ! allocate sds_vmr_inst to be used later. 
      ! for sufficient number of variables, make at least 7 more spaces available for 
      ! grid definition.
      allocate(sds_vmr_inst(nregions,nvmr,vmr_nvars_max+7))
      allocate(CFVMR_hdf_inst(nregions,nvmr))
      allocate(CFVMR_hdf_average(nregions,nvmr))
      allocate(sds_vmr_average(nregions,nvmr,vmr_nvars_max+7))
      
      do ivmr = 1, nvmr
        ! number in rc keys:
        write (nr,'(i3.3)') ivmr
        ! apply ?
        call ReadRc( rcF, 'output.cf.vmr.average.'//nr//'.apply', vmr_apply_average(ivmr), status )
        IF_NOTOK_RETURN(status=1)
        ! proceed ? for now use the data read in for instantaneous output...
        ! This has not been coded yet... see vmr
        if ( vmr_apply_average(ivmr) ) then
          ! first part of filename:
          ! call ReadRc( rcF, 'output.cf.vmr.average.'//nr//'.fname', vmr_fname_average(ivmr), status )
          ! IF_NOTOK_RETURN(status=1)
          ! time resolution:
          ! call ReadRc( rcF, 'output.cf.vmr.average.'//nr//'.dhour', vmr_dhour_average(ivmr), status )
          ! IF_NOTOK_RETURN(status=1)
          ! tracers to be written:
          ! call ReadRc( rcF, 'output.cf.vmr.average.'//nr//'.tracers', vmr_tracer_names_average(ivmr), status )
          ! IF_NOTOK_RETURN(status=1)
        end if  ! apply ?
      end do  ! vmr average numbers
    end if  ! nvmr > 0
    
    ! deposition fluxes: (not yet applicable... Please see 'vmr'-type of output)
    deps_fname = 'depositions'
    call ReadRc( rcF, 'output.cf.depositions.apply', deps_apply, status )
    IF_NOTOK_RETURN(status=1)
    if ( deps_apply ) then
      call ReadRc( rcF, 'output.cf.depositions.dhour', deps_dhour, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.depositions.tracers', deps_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! deposition velocities: (not yet applicable... Please see 'vmr'-type of output)
    depv_fname = 'depvels'
    call ReadRc( rcF, 'output.cf.depvels.apply', depv_apply, status )
    IF_NOTOK_RETURN(status=1)
    if ( depv_apply ) then
      call ReadRc( rcF, 'output.cf.depvels.dhour', depv_dhour, status )
      IF_NOTOK_RETURN(status=1)
      call ReadRc( rcF, 'output.cf.depvels.tracers', depv_tracer_names, status )
      IF_NOTOK_RETURN(status=1)
    end if
    
   
    ! no files open yet
    curr_day = -1
    
    ! lowest time frequency is 1 hour
    dhour_min = 1
    
    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine Output_CF_Init
  
  

  
  
  
  subroutine Output_CF_Step( region, idate_f, status )

    ! --- in/out ---------------------------------
    
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Output_CF_Step'
    
    ! --- begin -------------------------------
    
    call goLabel(rname)

    !
    ! close files if necessary
    !
    
    if ( all(curr_day(region,:) > 0) .and. any(idate_f(1:3) /= curr_day(region,:)) ) then
    
      ! write for 24:00
      call CF_Files_Write2( region,  idate_f, status )
      IF_NOTOK_RETURN(status=1)
      
      ! close all:
      call CF_Files_Close( region, status )
      IF_NOTOK_RETURN(status=1)
      
      ! no files open ...
      curr_day(region,:) = -1
  
    end if

    !
    ! open files if necessary
    !
    
    if ( any(curr_day(region,:) < 0) ) then
    
      ! open files:
      call CF_Files_Open( region, idate_f, status )
      IF_NOTOK_RETURN(status=1)
      
      ! store date:
      curr_day(region,:) = idate_f(1:3)
      
    end if

    !
    ! write files
    !
    
    call CF_Files_Write( region, idate_f, status )
    IF_NOTOK_RETURN(status=1)

    ! special files-write:

    ! if ( any(idate_f(4:6) > 0) ) then
    !   call CF_Files_Write2( region, idate_f, status )
    !   IF_NOTOK_RETURN(status=1)
    ! end if
    
    !
    ! done
    !
    
    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine Output_CF_Step
  
  
  ! ***
  
  
  subroutine Output_CF_Done( status )
  
    ! --- in/out ---------------------------------
    
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Output_CF_Done'
    
    ! --- local -------------------------------

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

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

    end do

    ! clear:
    if ( nvmr > 0 ) then
      deallocate( vmr_apply        )
      deallocate( vmr_apply_average)
      deallocate( vmr_fname        )
      deallocate( vmr_nvars        )
      deallocate( vmr_dhour        )
      deallocate( vmr_tracer_names )
      deallocate( CFVMR           )
      deallocate( CFVMR_hdf_inst)
      deallocate( CFVMR_hdf_average )
      deallocate(sds_vmr_inst)
      deallocate(sds_vmr_average)
    end if
      
    ! ok
    status = 0
    
  end subroutine Output_CF_Done
  
  
  ! ********************************************************************
  ! ***
  ! *** open/write/close retro files
  ! ***
  ! ********************************************************************
  

  subroutine CF_Files_Open( region, idate_f, status )
  
    use global_data, only : outdir

    ! --- in/out ---------------------------------
    
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/CF_Files_Open'
    
    ! --- local -------------------------------
    
    integer            ::  ivmr
    
    ! --- begin -------------------------------
    
    ! grid definition:
    if ( griddef_apply ) then
      call CFGridDef_Init( CFGridDef(region), GridDef_hdf(region), outdir, fname_model, fname_expid, region, status )
      IF_NOTOK_RETURN(status=1)
    end if

    ! dynamics:
!    if ( tp_apply ) then
!      call CFTP_Init     ( CFTP(region)     , outdir, fname_model, fname_expid, region, idate_f, status )
!      IF_NOTOK_RETURN(status=1)
!    end if

    ! tracer concentrations: inst/ mean
    do ivmr = 1, nvmr
      if ( vmr_apply(ivmr) ) then
      call CFVMR_Init(CFVMR_hdf_inst(region,ivmr), CFVMR(region,ivmr), outdir, fname_model, fname_expid, &
                            vmr_fname(ivmr), region, idate_f, &
                            vmr_dhour(ivmr), vmr_tracer_names(ivmr), & 
			    ivmr, status )
      IF_NOTOK_RETURN(status=1)
      endif
    end do

    ! deposition fluxes:
!    if ( deps_apply ) then
!      call CFDEPS_Init( CFDEPS(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 CFDEPV_Init( CFDEPV(region), outdir, fname_model, fname_expid, &
!                            depv_fname, region, idate_f, &
!                            depv_dhour, depv_tracer_names, status )
!      IF_NOTOK_RETURN(status=1)
!    end if


    ! ok
    status = 0
    
  end subroutine CF_Files_Open
  
  
  ! ***
  
  
  subroutine CF_Files_Write( region, idate_f, status )
  
    ! --- in/out ---------------------------------
    
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/CF_Files_Write'
    
    ! --- local -------------------------------
    
    integer            ::  ivmr
    
    ! --- begin -------------------------------
    
    ! grid definition:
    if ( griddef_apply ) then
      call CFGridDef_Write( CFGridDef(region), GridDef_hdf(region), region, status )
      IF_NOTOK_RETURN(status=1)
    end if
    
    ! dynamics:
!    if ( tp_apply ) then
!      call CFTP_Write( CFTP(region), region, idate_f, status )
!      IF_NOTOK_RETURN(status=1)
!    end if
    
    ! instantaneous tracer fields:
    do ivmr = 1, nvmr
      if ( .not. vmr_apply(ivmr) ) cycle
      call CFVMR_Write( CFVMR(region,ivmr), region, ivmr,idate_f, status )
      IF_NOTOK_RETURN(status=1)
    end do

    ! average tracer fields:
    ! do ivmr = 1, nvmr
    !   if ( .not. vmr_apply_average(ivmr) ) cycle
    !   call CFVMR_average( CFVMR(region,ivmr), region, idate_f, status )
    !   IF_NOTOK_RETURN(status=1)
    ! end do

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



    ! ok
    status = 0
    
  end subroutine CF_Files_Write
  
  
  ! ***
  
  
  ! write at end of time interval
  
  subroutine CF_Files_Write2( region, idate_f, status )
  
    ! --- in/out ---------------------------------
    
    integer, intent(in)              ::  region
    integer, intent(in)              ::  idate_f(6)
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/CF_Files_Write2'
    
    ! --- begin -------------------------------
    
    ! deposition fluxes:
!    if ( deps_apply ) then
!      call CFDEPS_Write( CFDEPS(region), region, idate_f, status )
!      IF_NOTOK_RETURN(status=1)
!    end if

    ! ok
    status = 0
    
  end subroutine CF_Files_Write2
  
  
  ! ***
  
  
  subroutine CF_Files_Close( region, status )
  
    ! --- in/out ---------------------------------
    
    integer, intent(in)              ::  region
    integer, intent(out)             ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/CF_Files_Close'
    
    ! --- local -------------------------------
    
    integer            ::  ivmr
    
    ! --- begin -------------------------------
    
    if ( griddef_apply ) then
      call CFGridDef_Done( CFGridDef(region),GridDef_hdf(region), status )
      IF_NOTOK_RETURN(status=1)
    end if

!    if ( tp_apply ) then
!      call CFTP_Done     ( CFTP(region)     , status )
!      IF_NOTOK_RETURN(status=1)
!    end if

    do ivmr = 1, nvmr
      if ( .not. vmr_apply(ivmr) ) cycle
      call CFVMR_Done( CFVMR(region,ivmr),CFVMR_hdf_inst(region,ivmr),region,ivmr, status )
      IF_NOTOK_RETURN(status=1)
    end do

!    if ( deps_apply ) then
!      call CFDEPS_Done( CFDEPS(region), status )
!      IF_NOTOK_RETURN(status=1)
!    end if

!    if ( depv_apply ) then
!      call CFDEPV_Done( CFDEPV(region), status )
!      IF_NOTOK_RETURN(status=1)
!    end if


    ! ok
    status = 0
    
  end subroutine CF_Files_Close
  
  
  
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE 1: Model horizontal grid definition 
  ! (longitude, latitude, size of gridbox [m2] ). 
  ! native vertical grid definition (hybrid level coefficients) and the 
  ! formula used to calculate pressure.
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    

  subroutine CFGridDef_Init( CF, CF_hdf, fdir, model, expid, region, status )
  
    use partools, only : myid, root, MPI_INFO_NULL
    use partools, only : localComm
    use Meteo   , only : lli, levi
    use file_hdf, only : THdfFile, Init,  setdim,WriteAttribute, WriteData, TSds, Done
    ! --- in/out -------------------------------------
    
    type(TCF_File_GridDef), intent(out)    ::  CF
    type(THdfFile),intent(out)             ::  CF_hdf
    character(len=*), intent(in)           ::  fdir
    character(len=*), intent(in)           ::  model
    character(len=*), intent(in)           ::  expid
    integer, intent(in)                    ::  region
    integer, intent(out)                   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFGridDef_Init'
    
    ! --- local ------------------------------------
    
!    type(THdfFile)        ::  CF_hdf
    character(len=256)    ::  fname
    integer               ::  j,varid,ivar
    character(len=256)    ::  xname

    ! --- begin -------------------------------------
    
    if (myid /= root) return 

    call goLabel(rname)
    
    ! o open file
    ! write filename
    write (fname,'(a,"/",a,a,"_",a,"_",a,".hdf")') &
             trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), 'griddef'
    call Init(CF_hdf,fname,'create', status)
    IF_NOTOK_RETURN(status=1)
    
    ! o global attributes

    call WriteAttribute(CF_hdf,'title','GridDef',status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(CF_hdf,'dataset_author', trim(dataset_author),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(CF_hdf,'dataset_institution', trim(dataset_institution),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(CF_hdf,'dataset_version', trim(dataset_version),status)
    IF_NOTOK_RETURN(status=1)

    ! o define dimensions
    
      CF%dimid_lon = lli(region)%nlon
      CF%dimid_lat = lli(region)%nlat
      CF%dimid_lev = levi%nlev
      CF%dimid_levi= levi%nlev+1


    ! o define variables

    ivar=1  ! longitude
    call init(sds_grid(ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lon /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_grid(ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units),(/(lli(region)%lon_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )
    IF_NOTOK_RETURN(status=1)

    ivar = 2 ! latitude
    call init(sds_grid(ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lat /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_grid(ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units),(/(lli(region)%lat_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )


    ivar = 6 ! gridbox area
    call init(sds_grid(ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lon ,CF%dimid_lat /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_grid(ivar), 0, trim(CF_grid_var(1)%param_name), trim(CF_grid_var(1)%units),(/(lli(region)%lon_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_grid(ivar), 1, trim(CF_grid_var(2)%param_name), trim(CF_grid_var(2)%units),(/(lli(region)%lat_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )


    do ivar = 7,8 ! a,b mid-level
    call init(sds_grid(ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lev /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_grid(ivar), 0, trim(CF_grid_var(3)%param_name), trim(CF_grid_var(3)%units),(/(j,j=1,levi%nlev)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )
    enddo
    
    do ivar = 9,10 ! a,b boundary values (at cell interfaces)
    call init(sds_grid(ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_levi /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_grid(ivar), 0, 'levi', trim(CF_grid_var(3)%units),(/(j,j=1,levi%nlev+1)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_grid(ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )
    enddo

 
    !  o end defintion mode

    
    !  o
    
    ! no records written yet
    CF%trec = 0
    
    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine CFGridDef_Init
  
  
  ! ***
  
    
  subroutine CFGridDef_Write( CF, CF_hdf, region, status )
  
    use GO      , only : TDate, NewDate, rTotal, operator(-)
    use Grid    , only : AreaOper
    use partools, only : myid, root
    use Meteo   , only : lli, levi
    use Meteo   , only : sp_dat
    use file_hdf,    only : THdfFile, WriteData
    ! --- in/out -------------------------------------
    
    type(TCF_File_GridDef), intent(inout)  ::  CF
    type(THdfFile), intent(inout)          ::  CF_hdf
    integer, intent(in)                      ::  region
    integer, intent(out)                     ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFGridDef_Write'
    
    ! --- local ------------------------------------
    
    integer             ::  imr, jmr, lmr,ivar
    real, allocatable   ::  ll(:,:)
    !type(TDate)         ::  t, t0
    real                ::  time
    
    ! --- begin -------------------------------------
    
    call goLabel(rname)
    
    ! grid size
    imr = lli(region)%nlon
    jmr = lli(region)%nlat
    lmr = levi%nlev
    
    ! next time record:
    CF%trec = CF%trec + 1
    
    !! time since 2000-1-1 00:00
    !t0 = NewDate( time6=time_reftime6 )
    !t = NewDate( time6=idate_f )
    !time = rTotal( t - t0, 'day' )
    
    
    ! root only:
    if ( myid == root ) then
    
      ! lat/lon field:
      allocate( ll(imr,jmr) )

      ! o write data
   
      if ( CF%trec == 1 ) then
      ivar=1
      call writedata(sds_grid(ivar),lli(region)%lon_deg,status)
      IF_NOTOK_RETURN(status=1)
      ivar=2
      call writedata(sds_grid(ivar),lli(region)%lat_deg,status)
      IF_NOTOK_RETURN(status=1)
 
      ll = 1.0
      call AreaOper( lli(region), ll, '*', 'm2', status )
      IF_NOTOK_RETURN(status=1)

      ivar=6
      call writedata(sds_grid(ivar),ll,status)
      IF_NOTOK_RETURN(status=1)
 
      ivar=7
      call writedata(sds_grid(ivar), levi%fa ,status)
      IF_NOTOK_RETURN(status=1)

      ivar=8
      call writedata(sds_grid(ivar),levi%fb ,status)
      IF_NOTOK_RETURN(status=1)

      ivar=9
      call writedata(sds_grid(ivar),levi%a(0:levi%nlev),status)
      IF_NOTOK_RETURN(status=1)

      ivar=10
      call writedata(sds_grid(ivar),levi%b(0:levi%nlev),status)
      IF_NOTOK_RETURN(status=1)

     end if


      ! clear
      deallocate( ll )
    
    end if   ! root
    
    ! end independend data mode:
    !status = pnf90_end_indep_data( CF%ncid )
    !IF_PNF90_NOTOK_RETURN(status=1)

    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine CFGridDef_Write
  

  ! ***
  
    
  subroutine CFGridDef_Done( CF, CF_hdf, status )
  
    use partools, only : myid, root
    use Meteo   , only : lli, levi
    use file_hdf, only : THdfFile,Done
  
    ! --- in/out -------------------------------------
    
    type(TCF_File_GridDef), intent(inout)  ::  CF
    type(THdfFile), intent(inout)          ::  CF_hdf
    integer, intent(out)                   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFGridDef_Done'
    
    ! --- local ------------------------------------
    integer           ::  ivar
    ! --- begin -------------------------------------
    
    call goLabel(rname)
    if (myid == root) then
    ! close file
    
    call Done(sds_grid(1),status)
    call Done(sds_grid(2),status)
    call Done(sds_grid(6),status)
    call Done(sds_grid(7),status)
    call Done(sds_grid(8),status)
    call Done(sds_grid(9),status)
    call Done(sds_grid(10),status)

    call Done(CF_hdf, status)
    IF_NOTOK_RETURN(status=1)

    end if   ! myid == root
    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine CFGridDef_Done


  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE2: 3D field of monthly Model pressure [Pa] and temperature [K].
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    

  subroutine CFTP_Init( CF, fdir, model, expid, region, idate_f, status )
  
    use partools, only : myid, root, MPI_INFO_NULL
    use partools, only : localComm
    use Meteo   , only : lli, levi
    use Meteo   , only : Set
    use Meteo   , only : sp_dat, oro_dat, temper_dat, humid_dat, pu_dat, pv_dat, mfw_dat, gph_dat, t2m_dat
  
    ! --- in/out -------------------------------------
    
    type(TCF_File_TP), intent(out)       ::  CF
    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(out)                   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFTP_Init'
    
    ! --- local ------------------------------------
    
    character(len=256)    ::  fname
    integer               ::  varid
    
    ! --- begin -------------------------------------
    
    call goLabel(rname)
    

    ! ensure that required meteo is loaded:
    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

    ! o open file
    
    ! write filename
!    write (fname,'(a,"/",a,a,"_",a,"_",a,"_",i4.4,"_",i2.2,"_",i2.2,".hdf")') &
!             trim(fdir), trim(model), trim(fname_grid(region)), trim(expid), 'TP', idate_f(1:3)
    
    ! open:
!    status = pnf90_create( localComm, fname, PNF90_CLOBBER, MPI_INFO_NULL, CF%ncid )
!    IF_PNF90_NOTOK_RETURN(status=1)
    
    ! o global attributes
    
    
    ! o define dimensions


    ! o define variables


    !  o end defintion mode

    !  o
    
    ! no records written yet
    CF%trec = 0
    
    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine CFTP_Init
  
  
  ! ***
  
    
  subroutine CFTP_Write( CF, region, idate_f, status )

    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 Meteo      , only : lli, levi
    use Meteo      , only : sp_dat, temper_dat, humid_dat, pu_dat, pv_dat, mfw_dat, gph_dat, oro_dat, t2m_dat
    use Meteo      , only : m_dat
    use global_data, only : mass_dat
  
    ! --- in/out -------------------------------------
    
    type(TCF_File_TP), intent(inout)     ::  CF
    integer, intent(in)                    ::  region
    integer, intent(in)                    ::  idate_f(6)
    integer, intent(out)                   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFTP_Write'
    
    ! --- local ------------------------------------
    
    integer             ::  i, j, l
    integer             ::  imr, jmr, lmr
    real, allocatable   ::  lev(:)
    type(TDate)         ::  t, t0
    real                ::  time
    real, allocatable   ::  field3d(:,:,:)
    real, allocatable   ::  p_hlev(:)
    
    ! --- begin -------------------------------------
    
    call goLabel(rname)
    
    ! grid size
    imr = lli(region)%nlon
    jmr = lli(region)%nlat
    lmr = levi%nlev
    
    ! next time record:
    CF%trec = CF%trec + 1
    
    ! time since reftime:
    t0 = NewDate( time6=time_reftime6 )
    t = NewDate( time6=idate_f )
    time = rTotal( t - t0, 'day' )
    
    ! start independend data mode:
      ! o write data
      
    ! end independend data mode:
    
    call goLabel()

    ! ok
    status = 0

  end subroutine CFTP_Write
  

  ! ***
  
    
  subroutine CFTP_Done( CF, status )
  
    use partools, only : myid, root
    use Meteo   , only : lli, levi
  
    ! --- in/out -------------------------------------
    
    type(TCF_File_TP), intent(inout)  ::  CF
    integer, intent(out)                     ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFTP_Done'
    
    ! --- local ------------------------------------
    
    ! --- begin -------------------------------------
    
    call goLabel(rname)
    
    ! close file
    
    call goLabel()
      
    ! ok
    status = 0
    
  end subroutine CFTP_Done


  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  ! FILE3: 3D fields for O3, CO, CH4, ...  volume mixing ratios
  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
  subroutine CFVMR_Init(CF_hdf, CF, fdir, model, expid, filetype, region, &
                              idate_f, dhour, tracer_names, ivmr, status )

    use Binas      , only : xmair
    use GO        , only : goReadFromLine, goUpCase
    use chem_param , only : ntrace, names, ra
    use partools   , only : myid, root, MPI_INFO_NULL
    use partools   , only : localComm
    use Meteo      , only : lli, levi, sp_dat, cp_dat
    use file_hdf,    only : THdfFile, Init,  setdim,WriteAttribute, WriteData, TSds, Done,SD_UNLIMITED

    ! --- in/out -------------------------------------
    
    type(THdfFile),intent(out)             ::  CF_hdf
    type(TCF_File_VMR), intent(out)        ::  CF
    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)                    ::  ivmr
    integer, intent(in)                    ::  idate_f(6)
    integer, intent(in)                    ::  dhour
    character(len=*), intent(in)           ::  tracer_names
    integer, intent(out)                   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFVMR_Init'
    
    ! --- local ------------------------------------
    
    character(len=256)    ::  fname
    integer               ::  varid
    
    character(len=256)    ::  trnames
    character(len=8)      ::  trname, tmname
    integer               ::  k, j,itr,ivar,ivar_tracer,ivar_vmr
    
    character(len=32)     ::  varname, varname_conc, 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
    CF%dhour        = dhour
    CF%tracer_names = tracer_names
    

    !
    ! Find trace index for requested tracers.
    !

    write (gol,'("selected tracers for VMR output:")'); call goPr
    CF%ntr = 0
    trnames = tracer_names
    do
      ! empty ?
      if ( len_trim(trnames) == 0 ) exit
      ! next number:
      if ( CF%ntr == ntrace ) then
        write (gol,'("number of elements in tracer names list exceeds ntrace=",i6)') ntrace; call goErr
        TRACEBACK; status=1; return
      end if
      CF%ntr = CF%ntr + 1
      ! extract leading name:
      call goReadFromLine( trnames, trname, status, sep=' ' )
      IF_NOTOK_RETURN(status=1)
      ! store cf name:
      CF%name_tr(CF%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
      ! NOy is a special ...
      select case ( tmname )
        case ( 'NOy' )
          ! defined as ntrace+1
          CF%itr(CF%ntr) = iNOy
          write (gol,'("    *  ",a10)') trim(trname); call goPr
        case default
          ! loop over all names:
          CF%itr(CF%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
              CF%itr(CF%ntr) = itr
              exit
            end if
          end do
      end select
      ! not found ?
      if ( CF%itr(CF%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)') CF%ntr; call goPr
        write (gol,'("  cf    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 ( CF%ntr < 1 ) then
      write (gol,'("no tracers extracted from list :",a)') tracer_names; call goErr
      TRACEBACK; status=1; return
    end if

    !
    ! Only root should initialize files
    !
    if (myid /= root) return 

    ! o open file for writing instantaneous data:
    ! write filename

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

    ! open:
    call Init(CF_hdf,trim(fname),'create', status)
    IF_NOTOK_RETURN(status=1)

    ! o write global attributes

    call WriteAttribute(CF_hdf,'title','instantaneous volume mixing ratios',status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(CF_hdf,'dataset_author', trim(dataset_author),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(CF_hdf,'dataset_institution', trim(dataset_institution),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(CF_hdf,'dataset_version', trim(dataset_version),status)
    IF_NOTOK_RETURN(status=1)

    ! o define dimensions

    CF%dimid_lon = lli(region)%nlon
    CF%dimid_lat = lli(region)%nlat
    CF%dimid_lev = levi%nlev
    !This doesn't work. I don't know what should be set here...:      CF%dimid_time='UNLIMITED'
    CF%dimid_time=2
    CF%dimid_datelen=6
     
    ! o define variables, using the CF_grid_var and CF_vmr_var -stuff defined earlier

    ivar=1  ! longitude
    call init(sds_vmr_inst(region,ivmr,ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lon /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units),(/(lli(region)%lon_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )
    IF_NOTOK_RETURN(status=1)

    ivar = 2 ! latitude
    call init(sds_vmr_inst(region,ivmr,ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lat /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units),(/(lli(region)%lat_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )

    ivar = 3 ! level
    call init(sds_vmr_inst(region,ivmr,ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_lev /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units),(/(j,j=1,levi%nlev)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )

    ivar = 4 ! time
    call init(sds_vmr_inst(region,ivmr,ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ SD_UNLIMITED /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
!    call SetDim( sds_vmr_inst(region,ivmr,ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units), (/ /) ,status)
!    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )
    IF_NOTOK_RETURN(status=1)

    ivar = 5 ! date
    call init(sds_vmr_inst(region,ivmr,ivar), CF_hdf, trim(CF_grid_var(ivar)%param_name), (/ CF%dimid_datelen, SD_UNLIMITED /), 'integer(4)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 0, trim(CF_grid_var(ivar)%param_name), trim(CF_grid_var(ivar)%units),(/1,2,3,4,5,6/) ,status)
    IF_NOTOK_RETURN(status=1)
!    call SetDim( sds_vmr_inst(region,ivmr,ivar), 1, trim(CF_grid_var(4)%param_name), trim(CF_grid_var(4)%units),(/0/) ,status)
!    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'standard_name'  , trim(CF_grid_var(ivar)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'long_name'      ,trim(CF_grid_var(ivar)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'units'          , trim(CF_grid_var(ivar)%units),status )
    IF_NOTOK_RETURN(status=1)

    ivar_vmr = 6  ! surface pressure 
    call init(sds_vmr_inst(region,ivmr,ivar_vmr), CF_hdf, trim(CF_TP_var(1)%param_name), (/CF%dimid_lon,CF%dimid_lat, SD_UNLIMITED /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar_vmr), 0, trim(CF_grid_var(1)%param_name), trim(CF_grid_var(1)%units),(/(lli(region)%lon_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar_vmr), 1, trim(CF_grid_var(2)%param_name), trim(CF_grid_var(2)%units),(/(lli(region)%lat_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar_vmr), 'standard_name'  , trim(CF_TP_var(1)%stand_name) ,status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar_vmr), 'long_name'      , trim(CF_TP_var(1)%long_name) ,status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar_vmr), 'units'          , trim(CF_TP_var(1)%param_name),status )
    IF_NOTOK_RETURN(status=1)


    nvar_dims=6

    ! loop over tracer to be written:
    do k = 1, CF%ntr
    ivar_tracer=CF%itr(k)

    ivar = nvar_dims+k

    call init(sds_vmr_inst(region,ivmr,ivar), CF_hdf, trim(CF_vmr_var(ivar_tracer)%param_name), (/CF%dimid_lon,CF%dimid_lat, CF%dimid_lev, SD_UNLIMITED /), 'real(8)',status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 0, trim(CF_grid_var(1)%param_name), trim(CF_grid_var(1)%units),(/(lli(region)%lon_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 1, trim(CF_grid_var(2)%param_name), trim(CF_grid_var(2)%units),(/(lli(region)%lat_deg)/) ,status)
    IF_NOTOK_RETURN(status=1)
    call SetDim( sds_vmr_inst(region,ivmr,ivar), 2, trim(CF_grid_var(3)%param_name), trim(CF_grid_var(3)%units),(/(j,j=1,levi%nlev)/) ,status)
    IF_NOTOK_RETURN(status=1)
!    call SetDim( sds_vmr_inst(region,ivmr,ivar), 3, trim(CF_grid_var(4)%param_name), trim(CF_grid_var(4)%units),(/0/) ,status)
!    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'standard_name'  , trim(CF_vmr_var(ivar_tracer)%stand_name),status)
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'long_name'      ,trim(CF_vmr_var(ivar_tracer)%long_name),status )
    IF_NOTOK_RETURN(status=1)
    call WriteAttribute(sds_vmr_inst(region,ivmr,ivar), 'units'          , trim(CF_vmr_var(ivar_tracer)%units),status )
    IF_NOTOK_RETURN(status=1)

      
      ! store varid (?)
      CF%varid_tr(k) = varid
      
    end do

    !  o end defintion mode

    ! status = pnf90_enddef( CF%ncid )
    ! IF_PNF90_NOTOK_RETURN(status=1)
    
    !  o
    
    ! no records written yet
    CF%trec = 0
    
    call goLabel()
    
    ! ok
    status = 0
    
  end subroutine CFVMR_Init
  
  
  ! ***
  
    
  subroutine CFVMR_Write( CF, region, ivmr, idate_f, status )
  
    use GO         , only : TDate, NewDate, rTotal, operator(-)
    use chem_param , only : ntrace, ntracet, fscale,names,ntrace_chem
    use partools   , only : myid, root
    use partools   , only : previous_par
    use partools   , only : tracer_active, tracer_loc, tracer_id
    use partools   , only : lmloc, offsetl
    use partools   , only : Par_Barrier, Par_Scatter_Over_Levels, Par_Gather_Tracer_t
    use partools   , only : which_par, previous_par,par_gather_from_levels,par_gather_tracer_k
    use tracer_data, only : mass_dat, chem_dat
    use Meteo      , only : lli, levi
    use Meteo      , only : m_dat, sp_dat, cp_dat, iwc_dat, temper_dat
    use file_hdf   , only : THdfFile, WriteData, TSds, Done
    use dims       , only : jsr,jer,isr,ier
!vh    use emission_nox, only: flashrate
#ifdef MPI 
      use mpi_comm,  only : gather_tracer_t,gather_tracer_k
      use mpi_const, only : my_real, mpi_sum, com_trac,com_lev, ierr
#endif
   
    ! --- in/out -------------------------------------
    
    type(TCF_File_VMR), intent(inout)      ::  CF
    integer, intent(in)                    ::  region
    integer, intent(in)                    ::  ivmr
    integer, intent(in)                    ::  idate_f(6)
    integer, intent(out)                   ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFVMR_Write'
    
    ! --- local ------------------------------------
    
    integer             ::  imr, jmr, lmr
    real, allocatable   ::  lev(:)
    integer             ::  i,j,l, ls, le,nsend,lmm,n
    type(TDate)         ::  t, t0
    real(8),dimension(1)::  time
    integer             ::  k, itr, itr_loc,ivar
    integer             ::  k_comp, itr_comp, itr_comp_loc,lglob
    integer,dimension(20)::  my_nontr
    integer             ::  my_nnontr  
 
    real, dimension(:,:,:,:), pointer       ::  rm
    real, dimension(:,:,:,:), pointer       ::  rm_c
    real, dimension(:,:,:),   pointer       ::  m
    real, dimension(:,:,:,:), allocatable   ::  rmk  
    real, dimension(:,:,:,:),   allocatable   ::  x_t  ,rm_cg,rm_tg,x_chem

    
    real, allocatable   ::  NOy_k(:,:,:)
    real, allocatable   ::  field_t(:,:,:)
    real, allocatable   ::  field_k(:,:,:)
    real, allocatable   ::  flsum(:,:)
    real, allocatable   ::  flavg(:,:)

    
    ! --- begin -------------------------------------
    
    ! for multiple of dhour only ...
    if ( (modulo(idate_f(4),CF%dhour)/=0) .or. any(idate_f(5:6)/=0)  ) then
      status=0; return
    end if
    
    call goLabel(rname)
    
    ! grid size
    imr = lli(region)%nlon
    jmr = lli(region)%nlat
    lmr = levi%nlev
    
    ! next time record:
    CF%trec = CF%trec + 1
    
    ! time since 2000-1-1 00:00
    t0 = NewDate( time6=time_reftime6 )
    t = NewDate( time6=idate_f )
    time = rTotal( t - t0, 'day' )
    ! start independend data mode:
    
    ! root only:
    if ( myid == root ) then
      ! only once ...    
      if ( CF%trec == 1 ) then

         ivar=1
         call writedata(sds_vmr_inst(region,ivmr,ivar),lli(region)%lon_deg,status)
         IF_NOTOK_RETURN(status=1)

         ivar=2
         call writedata(sds_vmr_inst(region,ivmr,ivar),lli(region)%lat_deg,status)
         IF_NOTOK_RETURN(status=1)

         ! write level indices:
         allocate( lev(lmr) )
         do l = 1, lmr
            lev(l) = real(l)
         end do
         ivar=3
         call writedata(sds_vmr_inst(region,ivmr,ivar),lev,status)
         IF_NOTOK_RETURN(status=1)
         deallocate(lev)
        
      end if  ! first record

      ! time:
      ivar=4
      call WriteData( sds_vmr_inst(region,ivmr,ivar), real(time), status, start=(/CF%trec-1/)  )
!      call writedata( sds_vmr_inst(region,ivmr,ivar), time, status  )
      IF_NOTOK_RETURN(status=1)
      
      ! date:
      ivar=5
      call WriteData( sds_vmr_inst(region,ivmr,ivar), reshape((idate_f),(/6,1/)), status, start=(/0,CF%trec-1/)  )
      IF_NOTOK_RETURN(status=1)

      ! surface pressure
      ivar=6  ! Now ivar corresponds to sds_vmr_inst
      call WriteData(sds_vmr_inst(region,ivmr,ivar), sp_dat(region)%data(1:imr,1:jmr,1), status, start=(/0,0,CF%trec-1/))
      IF_NOTOK_RETURN(status=1)

    end if !Root dir
 
      which_par=previous_par(region)
      if ( which_par /= 'tracer' ) then
        write (gol,'("Wrong type of parallelization : ",a)') which_par; call goErr
        TRACEBACK; status=1; return
      endif

      rm_c => chem_dat(region)%rm_k
      m => m_dat(region)%data

      ! gather transported species

      allocate(x_t(-1:imr+2,-1:jmr+2,lmr, ntracet))
      if ( which_par == 'tracer' ) then 
        rm => mass_dat(region)%rm_t
        call gather_tracer_t(x_t,imr,jmr,lmr,2,2,0,ntracet,rm,.false.)
        nullify(rm)   
      else
        rm => mass_dat(region)%rm_k
        call gather_tracer_k(x_t,imr,jmr,lmr,2,2,0,ntracet,rm,.false.)
        nullify(rm) 
      end if

!
! gather non-transported species
! These are always parallelized over layers.
! First set tracers at correct location
!

! are there any non-transported tracers in the requested output list? 
! put these tracers in 'my_nontr'
! and evaluate the total number in 'my_nnontr'
      my_nnontr=0
      do k = 1,CF%ntr
         if (CF%itr(k) >= ntracet+1) then
           my_nnontr=my_nnontr+1
           my_nontr(my_nnontr)=CF%itr(k)
         endif 
      enddo

    if (my_nnontr>0) then
      allocate(x_chem(1:imr,1:jmr,lmr,my_nnontr))

#ifdef MPI   
      allocate(rm_cg (1:imr,1:jmr,lmr,my_nnontr))
      if (lmloc > 0) then
         lmm = offsetl
         do n=1,my_nnontr
            do l=1,lmloc
               lglob=l+offsetl  !offset is zero on tracer domain
               do j = 1,jmr
                 do i = 1,imr
                    rm_cg(i,j,lglob,n) = rm_c(i,j,l,my_nontr(n))
                 end do
               end do
            end do
         end do
      endif    
      nsend=imr*jmr*lmr*my_nnontr
      call mpi_allreduce( rm_cg, x_chem, nsend, &
                          my_real, mpi_sum, com_trac, ierr )

      deallocate(rm_cg)
#else
      do n=1,my_nnontr
         do l=1,lmr
            do j = 1,jmr
              do i = 1,imr
                 x_chem(i,j,l,n) = rm_c(i,j,l,my_nontr(n))
              end do
            end do
         end do
      end do
#endif
    endif   ! (my_nnontr >0)

    my_nnontr=0
    ! loop over all tracer to be written:
    do k = 1, CF%ntr
    
      ! global tracer index:
      itr = CF%itr(k)
	           
      ! sds-index:
      ivar= nvar_dims+k
  
      ! transported or chemistry only ?
      if ( (itr >= 1) .and. (itr <= ntracet) ) then
      
        if ( myid==root) then	  
          call WriteData( sds_vmr_inst(region,ivmr,ivar), &
	  reshape( x_t(1:imr,1:jmr,1:lmr,itr)/m(1:imr,1:jmr,1:lmr)*fscale(itr), (/imr,jmr,lmr/) ), & 
	  status, start=(/0,0,0,CF%trec-1/)  )
          IF_NOTOK_RETURN(status=1)

        endif 
      else if ( (itr >= ntracet+1) .and. (itr <= ntrace) ) then
      my_nnontr=my_nnontr+1      
        if ( myid==root) then	         
          call WriteData( sds_vmr_inst(region,ivmr,ivar), &
          reshape( x_chem(1:imr,1:jmr,1:lmr,my_nnontr)/m(1:imr,1:jmr,1:lmr)*fscale(itr), (/imr,jmr,lmr/) ), &
                   status, start=(/0,0,0,CF%trec-1/)  )
          IF_NOTOK_RETURN(status=1)
       endif 


      ! some exceptions...
      else if ( itr == iNOy ) then
      write(*,*)'ERROR: CHECK IMPLEMENTATION!!! for NOy'
      stop
        ! mole fraction of NOy = sum of mole fractions of NOy components
      
        ! storage for sum of NOy components (distributed over levels):
        ! allocate( NOy_k(imr,jmr,lmloc) )

        ! 3d fields with all levels or local levels only:
        ! allocate( field_t(imr,jmr,lmr  ) )
        ! allocate( field_k(imr,jmr,lmloc) )

        ! loop over transported NOy components:
        NOy_k = 0.0
        do k_comp = 1, nNOyt
          ! global tracer index:
          itr_comp = iNOyt(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
          ! how distributed ?
          select case ( previous_par(region) )
            ! distributed over tracers:
            case ( 'tracer' )
              ! fill 3D field with tracer or zeros:
              if ( tracer_active(itr_comp) ) then
                ! local tracer index:
                itr_comp_loc = tracer_loc(itr_comp)
                ! fill volume mixing ratio:
                field_t = mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_comp_loc) &
                           / m_dat(region)%data(1:imr,1:jmr,1:lmr) * fscale(itr_comp)
              else
                field_t = 0.0
              end if
              ! scatter from process with requested tracer over slabs of layers on all processors:
              call Par_Scatter_Over_Levels( field_t, tracer_id(itr_comp), field_k, status )
              IF_NOTOK_RETURN(status=1)
            ! distributed over layers:
            case ( 'levels' )
              ! copy into target array:
              field_k = mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr_comp) &
                         / m_dat(region)%data(1:imr,1:jmr,ls:le) * fscale(itr_comp)
            ! error ...
            case default
              write (gol,'("unsupported par for distributing NOy fields : ",a)') previous_par(region); call goErr
              TRACEBACK; status=1; return
          end select
          ! add contribution of this NOy component:
          NOy_k = NOy_k + field_k
        end do

       write (gol,'("Please implement treatment of variable : ",a)') 'NOy'; call goErr
        TRACEBACK; status=1; return

      else
      
        write (gol,'("strange tracer index requested for output : ",i6)') itr; call goErr
        TRACEBACK; status=1; return
        
      end if
    end do   ! tracer
    
   deallocate(x_t)
   if (my_nnontr>0) then
      deallocate(x_chem)
   endif
	 
    
    ! end
    call goLabel()

    ! ok
    status = 0
    
  end subroutine CFVMR_Write
  

  ! ***
  
    
  subroutine CFVMR_Done( CF, CF_hdf, region,ivmr,status )
  
    use partools, only : myid, root
    use Meteo   , only : lli, levi
    use file_hdf, only : THdfFile,Done
     
    ! --- in/out -------------------------------------
    
    type(TCF_File_VMR), intent(inout)        ::  CF
    type(THdfFile), intent(inout)            ::  CF_hdf
    integer, intent(in)                      ::  region
    integer, intent(in)                      ::  ivmr
    integer, intent(out)                     ::  status
    
    ! --- const --------------------------------------
    
    character(len=*), parameter ::  rname = mname//'/CFVMR_Done'
    
    ! --- local ------------------------------------

    integer                                  :: ivar,k

    ! --- begin -------------------------------------
    
    call goLabel(rname)
       if (myid == root) then
    ! close file
    

    do k = 1, nvar_dims
      call Done(sds_vmr_inst(region,ivmr,k),status)
      IF_NOTOK_RETURN(status=1)
    enddo

    do k = 1, CF%ntr
       call Done(sds_vmr_inst(region,ivmr,nvar_dims+k),status)
       IF_NOTOK_RETURN(status=1)
    enddo

    call Done(CF_hdf, status)
    IF_NOTOK_RETURN(status=1)

    end if   ! myid == root
    
        
    call goLabel()
      
    ! ok
    status = 0
    
  end subroutine CFVMR_Done
  
  
  end module user_output_cf