ecearth3/sources/tm5mp/proj/cbm4/user_output_aerocom.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_AEROCOM
!
! !DESCRIPTION: 
! 
!   This module provides objects and methods for the AEROCOM-2 related 
!   diagnostics package. It is strongly related to the setup of the 
!   optics module. 
!   A list of parameters are defined for which daily averages are output 
!   (so far), result is one file per day containing 50 2d-fields and 
!   4 3d-fields. 
! 
!   Life cycle(s): 
!   1) output_aerocom_init: 
!      - if( newsrun ): - allocation and initialisation of output container
!                         (mixf%..., drydepos/wetdepos/emission, wdep_out)
!                       - initialisation of optics module
!                          - lookuptable
!                          - wavelengths
!      - open output file (for today), initialise it and write attributes
!   2) output_aerocom_step:
!      - gather tracer
!      - add current state of model to containers 
!        (weighted by time since last call, --> period of validity)
!         - do fluxes, surface concentrations and loads
!         - do optics
!             - aerocom_aopio_init: initialise and fill input parameters       
!             - optics::optics_calculate_aop: 
!                 calculate optical properties (per wavelength)
!                 (take into account additional splitting, facilitated in 
!                  optics::optics_calculate_aop) 
!             - fill containers for optical properties
!             - free optics input parameters
!             - fill aerocom output containers 
!   3) output_aerocom_write: 
!      - collect_budgets: 
!         - use current budget containers (budemi, buddep, buddry)
!         - compare with saved values and infer by that the increment
!         - add to temporary output containers (wetdepos,drydepos,emission)
!      - write output containers to file 
!      - close file 
!   4) output_aerocom_done: 
!      - free allocated containers 
!        (mixf%..., drydepos/wetdepos/emission, wdep_out)
! 
!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_AEROCOM
  !
  ! !USES:
  !
  use go,           only : gol, goErr, goPr, goLabel
  use dims,         only : nregions
  use optics,       only : wavelendep
  use meteodata,    only : global_lli, levi
  use MDF
  use TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo, update_halo_iband

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: output_aerocom_init
  public :: output_aerocom_step
  public :: output_aerocom_write
  public :: output_aerocom_done
  public :: wdep_out
  character(len=20), public :: aerocom_exper ! AeroCom experiment name
  character(len=20), public :: aerocom_freq  ! AeroCom output frequency
  logical, parameter :: aai_output=.false.
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*),  parameter ::  mname = 'user_output_aerocom'

  ! Changed naming convections to AeroCom Control 2015 experiment
  character(len=20), parameter :: f_dataid='aerocom3', f_model='TM5'
  character(len=20), parameter :: f_dimmix='global',   f_dimstat='stations' 
  !
  ! !PRIVATE TYPES:
  !
  type metafields
     integer                       :: itm5
     character(len=16)             :: vname 
     character(len=64)             :: lname 
     character(len=10)             :: unit
     character(len=10)             :: positive
     character(len=130)            :: standard_name
  end type metafields

  type metafields2
     integer                       :: itm5
     character(len=17)             :: vname 
     character(len=64)             :: lname 
     character(len=10)             :: unit
     character(len=130)            :: standard_name
  end type metafields2

  type field3d 
     type( metafields )              :: mf
     real, dimension(:,:,:), pointer :: field
     integer                         :: varid
  end type field3d

  type field2d
     type( metafields )              :: mf
     real, dimension(:,:), pointer   :: field
     integer                         :: varid
  end type field2d

  type field1d
     type( metafields2 )             :: mf
     real, dimension(:,:), pointer   :: field
     integer                         :: varid
  end type field1d

  type field0d
     type( metafields2 )             :: mf
     real, dimension(:), pointer     :: field
     integer                         :: varid
  end type field0d

  type mixfile
     type( field3d ), dimension(:), pointer   :: f3d
     type( field2d ), dimension(:), pointer   :: f2d
     character(len=200)    :: fname
     integer               :: acct             ! accumulation time
     integer               :: funit            ! unit number 
     integer               :: nlon             ! x dimension of requested field
     integer               :: nlat             ! y dimension of requested field
     integer               :: nlev             ! z dimension of requested field
  end type mixfile
  type(mixfile), dimension(nregions), save  :: mixf  

  type stationfile
     type( field1d ), dimension(:), pointer   :: f1d
     type( field0d ), dimension(:), pointer   :: f0d
     character(len=200)    :: fname
     integer               :: acct             ! accumulation time
     integer               :: funit            ! unit number 
     integer               :: nstat            ! dimension of requested field
     integer               :: nlev             ! z dimension of requested field
  end type stationfile
  type(stationfile), dimension(nregions), save  :: statf  
  integer, dimension(5) :: statvarid

  type budgetstore
     real, dimension(:,:,:), allocatable :: f2dslast
     integer                             :: lasttime
  end type budgetstore

  type(budgetstore), dimension(nregions), save :: drydepos, wetdepos, emission

  ! wavelength information
  type(wavelendep), dimension(:), allocatable :: wdep_out 

  integer, parameter  :: statnamelen = 22
  integer             :: snamelendim
  type station
     integer           :: statnb                     ! station number
     character(len=statnamelen) :: statname          ! station name
     real              :: statlat                    ! geographic coordinate
     real              :: statlon                    ! geographic coordinate
     real              :: statalt                    ! geographic coordinate
  end type station
  type(station), dimension(:), allocatable :: statlist
  character(len=statnamelen), dimension(:), allocatable :: tmpstatnames

  type statintp
     integer,dimension(:),pointer :: ixr,  ixl  ! coordinates in current grid
     integer,dimension(:),pointer :: jyr,  jyl  ! coordinates in current grid
     integer,dimension(:),pointer :: lzr,  lzl  ! coordinates in current grid
     real,   dimension(:),pointer :: wixl, wixr ! interpolation coefficients
     real,   dimension(:),pointer :: wjyl, wjyr ! interpolation coefficients
     real,   dimension(:),pointer :: wlzl, wlzr ! interpolation coefficients
     real,   dimension(:),pointer :: tght       ! terrain height at station
  end type statintp
  type(statintp), dimension(nregions) :: sintp

  logical, dimension(:), allocatable :: stat_active    ! true if station is in processor domain
  integer                            :: stat_halo_loc  ! 1 if halo cells are needed for the local processor, else 0
  integer                            :: stat_halo_glob ! 1 if halo cells are needed on any processor, else 0
  ! the region dimension has already been removed for the above three variables

  integer                         ::  fid ! file id for IF_NOTOK_MDF macro
  integer                         ::  access_mode ! netcdf-4 parallel access mode for distributed data
  ! generate output if true:
  logical :: okdebug = .true.

  ! reference time:
  integer, parameter            ::  time_reftime6(6) = (/2001,01,01,00,00,00/)
  character(len=*), parameter   ::  time_units       = 'days since 2001-01-01 00:00' 

  ! sizes of arrays
  integer, save :: ntracer_3d, ntracer_2d
  integer, save :: ntracer_1dstat, ntracer_0dstat, nstations 
  integer, save :: nextra_1dstat

  ! index pointers for 1d/2d/3d fields in mixf
  integer, save :: temp, hus, airmass, pressure
  integer, save :: ec5503Daer, abs5503Daer, ec3503Daer, abs3503Daer
  integer, save :: ps          , precip      , sconcoa     , sconcbc     , sconcso4    , sconcdust   
  integer, save :: sconcss     , sconcno3    , loadoa      , loadbc      , loadso4     , loaddust    
  integer, save :: loadss      , loadno3     , emioa       , emibc       , emiso2      , emiso4      
  integer, save :: emidust     , emidms      , emiss       , dryso2      , drybc       , dryso4      
  integer, save :: drydust     , drydms      , dryss       , wetoa       , wetbc       , wetso2      
  integer, save :: wetso4      , wetdust     , wetdms      , wetss       , od550aer    , od550so4    
  integer, save :: od550bc     , od550oa     , od550ss     , od550dust   , od550no3    , od550aerh2o
  integer, save :: od550lt1aer , od550lt1dust, od550lt1ss  , abs550aer   , ec550aer    , asyaer
  integer, save :: ec550dryaer , abs550dryaer, asydryaer   , ec550drylt1aer, abs550drylt1aer
  integer, save :: od440aer    , od870aer    , sconcmsa    , dryoa       , sconcnh4
  integer, save :: abs440aer   , ec440dryaer , abs440dryaer
  integer, save :: abs870aer   , ec870dryaer , abs870dryaer
  integer, save :: od350aer    , abs350aer
  integer, save :: ec550dryaer_stat, abs550dryaer_stat, asydryaer_stat
  integer, save :: ec550drylt1aer_stat, abs550drylt1aer_stat
  integer, save :: ec440dryaer_stat, abs440dryaer_stat
  integer, save :: ec870dryaer_stat, abs870dryaer_stat
  integer, save :: bso4nus, bso4ais, bso4acs, bso4cos, bbcais , bbcacs 
  integer, save :: bbccos , bbcaii , bpomais, bpomacs, bpomcos, bpomaii
  integer, save :: bssacs , bsscos , bduacs , bducos , bduaci , bducoi, bno3_a 
  integer, save :: bnus_n , bais_n , bacs_n , bcos_n , baii_n , baci_n 
  integer, save :: bcoi_n , bh2onus, bh2oais, bh2oacs, bh2ocos
  !
  ! !REVISION HISTORY: 
  !   29 Nov 2010 - Achim Strunk   - Creation
  !   10 Oct 2014 - Twan van Noije - Adapted for TM5-mp
  !
  ! !REMARKS:
  ! (1) compiled only if with_m7 is used.
  !
  !EOP
  !------------------------------------------------------------------------

contains


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_AEROCOM_INIT
  !
  ! !DESCRIPTION: Initialise various parameters, eg., station lists,
  !               wavelengths for optics and output containers. 
  !               This routine should be called once per day, since fields
  !               are daily averages and one file per day is written.
  !               Additional tasks are done for newsrun==.true. . 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_aerocom_init(stat_output, status, iregion)
    !
    ! !USES:
    !
    !use MeteoData,     only : sp_dat, set
    use chem_param,    only : ntracet
    use chem_param,    only : inus_n, iso4nus, nmod
    use chem_param,    only : iais_n, iso4ais, ibcais, ipomais 
    use chem_param,    only : iacs_n, iso4acs, ibcacs, ipomacs, issacs, iduacs
    use chem_param,    only : icos_n, iso4cos, ibccos, ipomcos, isscos, iducos
    use chem_param,    only : iaii_n,          ibcaii, ipomaii
    use chem_param,    only : iaci_n, icoi_n,                           iduaci, iducoi
    use chem_param,    only : ino3_a, inh4, imsa
    use dims,          only : nregions, newsrun, idatee, idatei
    use dims,          only : region_name, parent, itau
    use dims,          only : xbeg, xend, ybeg, yend, dx, dy, xref, yref
    use dims,          only : zbeg, zend, dz, zref
    use global_data,   only : region_dat, outdir
    use datetime,      only : tau2date, date2tau
    use budget_global, only : nbud_vg
    use partools,      only : MPI_INFO_NULL, localComm
    use optics,        only : Optics_Init
    use datetime,      only : date2tau
    use dims,          only : mlen
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          :: status
    !
    ! !INPUT PARAMETERS:
    !
    logical, intent(in)           :: stat_output    ! true if stations output is requested
    integer, intent(in), optional :: iregion 
    !
    ! !REVISION HISTORY: 
    !   29 Nov 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer               :: imr, jmr, lmr, access_mode_sta
    integer               :: region, varid
    integer               :: io, istat, i, j, n, sc
    integer               :: i1, i2, j1, j2
    integer, dimension(6) :: idater
    character(len=10)     :: idates
    character(len=16)     :: lidates
    character(len=3)      :: cwavel
    real                  :: reftime
    integer               :: itaucur, itauref
    real                  :: dlat, dlon
    integer               :: mlength

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

    call goLabel(rname)

    ! access mode for distributed data (2d and 3d)
#ifdef MPI
#ifdef with_netcdf4_par       
    access_mode = MDF_COLLECTIVE
#else
    write(gol,'("AeroCom aerosol output requires netcdf4 with parallel access enabled")') ; call goErr
    TRACEBACK
    status=1; return
#endif
#else
    access_mode = MDF_INDEPENDENT
#endif

    ! for station
    access_mode_sta = MDF_INDEPENDENT

    !  initialize (only once)
    if( newsrun ) then 

       ! ensure that required meteo is loaded:
       ! call Set( sp_dat(region), status, used=.true. )

       ! set wavelength information 
       ! wl:    wavelength in microns
       ! split: whether to split into contributions from 
       !        M7 constituents (incl. water)
       allocate( wdep_out( 4 ) )
       wdep_out(1)%wl = 0.550 ; wdep_out(1)%split = .true.  ; wdep_out(1)%insitu = .true. 
       wdep_out(2)%wl = 0.440 ; wdep_out(2)%split = .false. ; wdep_out(2)%insitu = .true. 
       wdep_out(3)%wl = 0.870 ; wdep_out(3)%split = .false. ; wdep_out(3)%insitu = .true.
       wdep_out(4)%wl = 0.350 ; wdep_out(4)%split = .false. ; wdep_out(4)%insitu = .false.

       ! get the optics code prepared
       call Optics_Init(size(wdep_out), wdep_out, status )
       IF_NOTOK_RETURN(status=1)

       ! -----------------------
       ! parameters needed to reference the different 1d/2d/3d-fields
       ! (in order to avoid errors in referencing)
       ! finally, this list here simply determines the order in the output files
       ntracer_3d = 4
       if (aai_output) then
         ntracer_3d = ntracer_3d + 4
       endif
       ! -----------------------
       temp         =  1  ! 3d
       hus          =  2  ! 3d
       airmass      =  3  ! 3d
       pressure     =  4  ! 3d
       ec5503Daer   =  5  ! 3d
       abs5503Daer  =  6  ! 3d
       ec3503Daer   =  7  ! 3d
       abs3503Daer  =  8  ! 3d

       ntracer_2d   = 66
       ps           =  1  ! 2d
       precip       =  2  ! 2d
       sconcoa      =  3  ! 2d
       sconcbc      =  4  ! 2d
       sconcso4     =  5  ! 2d
       sconcdust    =  6  ! 2d
       sconcss      =  7  ! 2d
       sconcno3     =  8  ! 2d
       sconcnh4     =  9  ! 2d
       sconcmsa     = 10  ! 2d
       loadoa       = 11  ! 2d
       loadbc       = 12  ! 2d
       loadso4      = 13  ! 2d
       loaddust     = 14  ! 2d
       loadss       = 15  ! 2d
       loadno3      = 16  ! 2d
       emioa        = 17  ! 2d
       emibc        = 18  ! 2d
       emiso2       = 19  ! 2d
       emiso4       = 20  ! 2d
       emidust      = 21  ! 2d
       emidms       = 22  ! 2d
       emiss        = 23  ! 2d
       dryso2       = 24  ! 2d
       dryoa        = 25  ! 2d 
       drybc        = 26  ! 2d
       dryso4       = 27  ! 2d
       drydust      = 28  ! 2d
       drydms       = 29  ! 2d
       dryss        = 30  ! 2d
       wetoa        = 31  ! 2d
       wetbc        = 32  ! 2d
       wetso2       = 33  ! 2d
       wetso4       = 34  ! 2d
       wetdust      = 35  ! 2d
       wetdms       = 36  ! 2d
       wetss        = 37  ! 2d
       ! --- from here onwards keep consistent with order in optics (optics.F90)
       ! --- begin split order
       od550aer     = 38  ! 2d
       od550so4     = 39  ! 2d
       od550bc      = 40  ! 2d
       od550oa      = 41  ! 2d
       od550ss      = 42  ! 2d
       od550dust    = 43  ! 2d
       od550no3     = 44  ! 2d
       od550aerh2o  = 45  ! 2d
       od550lt1aer  = 46  ! 2d
       od550lt1dust = 47  ! 2d
       od550lt1ss   = 48  ! 2d
       ! --- end split order
       abs550aer    = 49  ! 2d
       asyaer       = 50  ! 2d
       ec550aer     = 51  ! 2d
       ! --- begin in-situ data order
       ec550dryaer  = 52  ! 2d
       abs550dryaer = 53  ! 2d
       asydryaer    = 54  ! 2d
       ec550drylt1aer  = 55  ! 2d
       abs550drylt1aer = 56  ! 2d
       ! --- end in-situ data order
       !
       od440aer     = 57  ! 2d
       abs440aer    = 58  ! 2d
       ec440dryaer  = 59  ! 2d
       abs440dryaer = 60  ! 2d
       !
       od870aer     = 61  ! 2d
       abs870aer    = 62  ! 2d
       ec870dryaer  = 63  ! 2d
       abs870dryaer = 64  ! 2d
       !
       od350aer     = 65  ! 2d
       abs350aer    = 66  ! 2d

       if (stat_output) then
          ! -----------------------
          ! set station information
          ! -----------------------
          ! this here is a kind of redundant coding, but it enables to easily check for correct 
          ! information (instead of assigning dedicated arrays (name,lat,lon etc.) side by side)

          ec550dryaer_stat=33
          abs550dryaer_stat=34
          asydryaer_stat=35
          ec550drylt1aer_stat=36
          abs550drylt1aer_stat=37
          ec440dryaer_stat=38
          abs440dryaer_stat=39
          ec870dryaer_stat=40
          abs870dryaer_stat=41
 
          nextra_1dstat=3 ! temperature, specific humidity and pressure at model levels
          ntracer_1dstat=41+nextra_1dstat
          ntracer_0dstat=ntracer_1dstat+2 ! plus surface pressure and precipitation

          ! Adapted station list for AeroCom Phase III in-situ projects
          nstations  = 90
          allocate( statlist( nstations ) )
          sc = 1
          statlist(sc) = station( sc,  'Alert                 ',  82.499,   -62.342,   210.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Anmyeon-do            ',  36.540,   126.330,    94.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Annaberg-Buchholz     ',  50.570,    13.000,   545.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Appalachian State U.  ',  36.213,   -81.692,  1076.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Aspvreten             ',  58.806,    17.388,    20.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Auchencorth Moss      ',  55.793,    -3.245,   260.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Backgarden            ',  23.540,   113.060,  -999.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Barrow                ',  71.323,  -156.611,    11.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'BEO Moussala          ',  42.179,    23.586,  2925.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Birkenes              ',  58.388,     8.252,   190.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Bondville             ',  40.050,   -88.367,   213.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Bukit Kototabang      ',  -0.202,   100.318,   864.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Bosel                 ',  52.998,     7.943,    40.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Cabauw                ',  51.971,     4.927,     1.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Cape Cod              ',  40.070,   -70.200,    43.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Cape Grim             ', -40.682,   144.688,    94.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Cape Point            ', -34.353,    18.490,   230.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Cape San Juan         ',  18.381,   -65.618,    65.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Chacaltaya            ', -16.200,   -68.100,  5320.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Danum Valley          ',   4.981,   117.844,   426.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Demokritos            ',  37.995,    23.816,   270.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'East Trout Lake       ',  54.350,  -104.983,   500.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Egbert                ',  44.230,   -79.783,   255.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'El Arenosillo         ',  37.104,    -6.734,    41.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'El Tololo             ', -30.173,   -70.799,  2220.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Elandsfontein         ', -25.533,    25.750,  1424.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Finokalia             ',  35.333,    25.670,   250.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Giordan Lighthouse    ',  36.073,    14.219,   167.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Gosan                 ',  33.280,   126.170,    72.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Graciosa              ',  39.080,   -28.030,    30.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Granada               ',  37.164,    -3.605,   680.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Hesselbach            ',  48.540,     8.400,    51.4 ); sc = sc + 1
          statlist(sc) = station( sc,  'Harwell               ',  51.573,    -1.317,   137.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Hohenpeissenberg      ',  47.802,    11.010,   985.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Hyytiala              ',  61.847,    24.295,   181.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Illmitz               ',  47.767,    16.767,   117.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Ispra                 ',  45.803,     8.627,   209.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Izana                 ',  28.309,   -16.499,  2373.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Jungfraujoch          ',  46.547,     7.985,  3580.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Kosetice              ',  49.583,    15.083,   535.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'K-Puszta              ',  46.967,    19.583,   125.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Leibzig               ',  51.353,    12.435,   126.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Leibzig-West          ',  51.318,    12.298,   122.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Lille Valby           ',  55.683,    12.117,    20.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'London - N. Kensington',  51.521,    -0.214,    27.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Lulin                 ',  23.470,   120.870,  2862.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Mace Head             ',  53.326,    -9.899,     5.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Madrid                ',  40.450,    -3.720,   669.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Manacapuro            ',  -3.210,   -60.590,    50.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Manaus                ',  -2.595,   -60.209,    45.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Mauna Loa             ',  19.536,  -155.576,  3397.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Melpitz               ',  51.530,    12.930,    86.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Montseny              ',  41.767,     2.350,   700.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Monte Cimone          ',  44.167,    10.683,  2165.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Mount Waliguan        ',  36.283,   100.900,  3810.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Nepal Climate Obs.    ',  27.958,    86.815,  5079.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Nainital              ',  29.360,    79.460,  1936.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Neuglobsov            ',  53.143,    13.033,    62.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Neumayer              ', -70.666,    -8.266,    42.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Niamey                ',  13.470,     2.170,   223.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Pallas                ',  67.974,    24.116,   560.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Payerne               ',  46.817,     6.950,   490.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Prague-Suchdol        ',  50.130,    14.380,   270.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Preila                ',  55.350,    21.067,     5.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Pt. Reyes             ',  38.091,  -122.957,     8.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Puy de Dome           ',  45.772,     2.966,  1465.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Resolute Bay          ',  74.717,   -94.983,    64.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Sable Island          ',  43.933,   -60.017,     2.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Samoa                 ', -14.232,  -170.563,    77.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Schauinsland          ',  47.900,     7.917,  1205.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Shang Dianzi          ',  40.390,   117.070,   294.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Shouxian              ',  32.558,   116.782,    23.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'SIRTA                 ',  48.709,     2.159,   160.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Sonnblick             ',  47.054,    12.959,  3106.0 ); sc = sc + 1 
          statlist(sc) = station( sc,  'South Pole            ', -89.997,   -24.800,  2841.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Southern Great Plains ',  36.617,   -97.500,   318.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Storm Peak            ',  40.455,  -106.744,  3220.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Summit                ',  72.580,   -38.480,  3238.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Tahkuse               ',  58.520,    24.940,    23.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Tiksi                 ',  71.586,   128.919,     8.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Trinidad Head         ',  41.054,  -124.151,   107.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Trollhaugen           ', -72.012,     2.535,  1553.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Uto                   ',  59.783,    21.383,     7.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Varrio                ',  67.767,    29.583,   400.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Vavihill              ',  56.017,    13.150,   172.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Vielsalm              ',  50.304,     6.001,   496.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Waldhof               ',  52.802,    10.759,    74.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Whistler Mountain     ',  50.059,  -122.958,  2182.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Zeppelin              ',  78.907,    11.889,   475.0 ); sc = sc + 1
          statlist(sc) = station( sc,  'Zugspitze             ',  47.417,    10.980,  2650.0 ); sc = sc + 1

          ! this copy is needed for the netcdf output of the 2d-character-array [statnamelen,nstations]
          allocate( tmpstatnames( nstations ) )
          tmpstatnames = statlist(:)%statname

          allocate(stat_active(nstations))
          stat_active(:)=.false.
          stat_halo_loc=0
       endif

    end if

    regionloop: do region = 1, nregions

       ! if region given, cycle if other region!
       if (present(iregion)) then
          if(iregion /= region) cycle regionloop
       endif

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

       if( nbud_vg /= lmr ) then
          write(gol,*)'output_aerocom_init: nbud_vg /= lmr'; call goErr
          write(gol,*)'output_aerocom_init: YOU MUST ADD THE PROJ/BUDGET10 TO SRC CODE !!!!!'; call goErr
          TRACEBACK
          status=1; return
       end if

       !  initialize the output:
       if( newsrun ) then 
 
          call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

          ! target structures for output
          allocate(mixf (region)%f3d(ntracer_3d))
          allocate(mixf (region)%f2d(ntracer_2d))

          ! allocate fields
          do i=1,ntracer_3d
             allocate( mixf(region)%f3d(i)%field(i1:i2,j1:j2,lmr) )
          end do
          do i=1,ntracer_2d
             allocate( mixf(region)%f2d(i)%field(i1:i2,j1:j2) )
          end do

          if (stat_output) then
             ! ---------
             ! find coordinates in grid and get interpolation coefficients
             ! ---------
             allocate( sintp(region)%ixl(nstations), sintp(region)%wixl(nstations), &
                       sintp(region)%ixr(nstations), sintp(region)%wixr(nstations), &
                       sintp(region)%jyl(nstations), sintp(region)%wjyl(nstations), &
                       sintp(region)%jyr(nstations), sintp(region)%wjyr(nstations), &
                       sintp(region)%lzl(nstations), sintp(region)%wlzl(nstations), &
                       sintp(region)%lzr(nstations), sintp(region)%wlzr(nstations), &
                       sintp(region)%tght(nstations) )
             call aerocom_get_statintp( region, status )
             IF_NOTOK_RETURN(status=1)
             allocate(statf(region)%f1d(ntracer_1dstat)) 
             allocate(statf(region)%f0d(ntracer_0dstat))
             do i=1,ntracer_1dstat
                allocate( statf(region)%f1d(i)%field(nstations,lmr) )
             end do
             do i=1,ntracer_0dstat
                allocate( statf(region)%f0d(i)%field(nstations) )
             end do 
          endif

          ! 3d data
          mixf(region)%f3d(temp        )%mf  = metafields( -1     , 'temp        ', 'Temperature                        ' , 'K         ', '', 'air_temperature'  )
          mixf(region)%f3d(hus         )%mf  = metafields( -1     , 'hus         ', 'Specific Humidity                  ' , '1         ', '', 'specific_humidity'  )
          mixf(region)%f3d(airmass     )%mf  = metafields( -1     , 'airmass     ', 'Air Mass                           ' , 'kg m-2    ', '', 'atmosphere_mass_of_air_per_unit_area'  )
          mixf(region)%f3d(pressure    )%mf  = metafields( -1     , 'pressure    ', 'Pressure                           ' , 'Pa        ', '', 'air_pressure'  )
          if (aai_output) then
            mixf(region)%f3d(ec5503Daer    )%mf  = metafields( -1   , 'ec5503Daer  ', 'Aerosol Extinction @550nm          ' , 'm-1       ', '',     'atmosphere_extinction_due_to_ambient_aerosol'  )
            mixf(region)%f3d(abs5503Daer   )%mf  = metafields( -1   , 'abs5503Daer ', 'Aerosol Absorption @550nm          ' , 'm-1       ', '',     'atmosphere_absorption_due_to_ambient_aerosol'  )
            mixf(region)%f3d(ec3503Daer    )%mf  = metafields( -1   , 'ec3503Daer  ', 'Aerosol Extinction @350nm          ' , 'm-1       ', '',     'atmosphere_extinction_due_to_ambient_aerosol'  )
            mixf(region)%f3d(abs3503Daer   )%mf  = metafields( -1   , 'abs3503Daer ', 'Aerosol Absorption @350nm          ' , 'm-1       ', '',     'atmosphere_absorption_due_to_ambient_aerosol'  )
          endif
          
          ! 2d data                    
          mixf(region)%f2d(ps          )%mf  = metafields( -1     , 'ps          ', 'Surface Air Pressure               ' , 'Pa        ', '',     'surface_air_pressure'  )
          mixf(region)%f2d(precip      )%mf  = metafields( -1     , 'precip      ', 'Precipitation                      ' , 'kg m-2 s-1', '',     'precipitation_flux'  )
          mixf(region)%f2d(sconcoa     )%mf  = metafields( -1     , 'sconcoa     ', 'Surface Concentration POM          ' , 'kg m-3    ', '',     'mass_concentration_of_particulate_organic_matter_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcbc     )%mf  = metafields( -1     , 'sconcbc     ', 'Surface Concentration BC           ' , 'kg m-3    ', '',     'mass_concentration_of_black_carbon_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcso4    )%mf  = metafields( -1     , 'sconcso4    ', 'Surface Concentration SO4          ' , 'kg m-3    ', '',     'mass_concentration_of_sulfate_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcdust   )%mf  = metafields( -1     , 'sconcdust   ', 'Surface Concentration Dust         ' , 'kg m-3    ', '',     'mass_concentration_of_dust_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcss     )%mf  = metafields( -1     , 'sconcss     ', 'Surface Concentration SS           ' , 'kg m-3    ', '',     'mass_concentration_of_seasalt_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcno3    )%mf  = metafields( -1     , 'sconcno3    ', 'Surface Concentration NO3          ' , 'kg m-3    ', '',     'mass_concentration_of_nitrate_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcnh4    )%mf  = metafields( -1     , 'sconcnh4    ', 'Surface Concentration NH4          ' , 'kg m-3    ', '',     'mass_concentration_of_ammonium_dry_aerosol_in_air'  )
          mixf(region)%f2d(sconcmsa    )%mf  = metafields( -1     , 'sconcmsa    ', 'Surface Concentration MSA          ' , 'kg m-3    ', '',     'mass_concentration_of_methane_sulfonic_acid_in_air'  )
          mixf(region)%f2d(loadoa      )%mf  = metafields( -1     , 'loadoa      ', 'Load of POM                        ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol'  )
          mixf(region)%f2d(loadbc      )%mf  = metafields( -1     , 'loadbc      ', 'Load of BC                         ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_black_carbon_dry_aerosol'  )
          mixf(region)%f2d(loadso4     )%mf  = metafields( -1     , 'loadso4     ', 'Load of SO4                        ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_sulfate_dry_aerosol'  )
          mixf(region)%f2d(loaddust    )%mf  = metafields( -1     , 'loaddust    ', 'Load of Dust                       ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_dust_dry_aerosol'  )
          mixf(region)%f2d(loadss      )%mf  = metafields( -1     , 'loadss      ', 'Load of SS                         ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_seasalt_dry_aerosol'  )
          mixf(region)%f2d(loadno3     )%mf  = metafields( -1     , 'loadno3     ', 'Load of NO3                        ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_nitrate_dry_aerosol'  )
          mixf(region)%f2d(emioa       )%mf  = metafields( -1     , 'emioa       ', 'Total Emission of POM              ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_due_to_emission'  )
          mixf(region)%f2d(emibc       )%mf  = metafields( -1     , 'emibc       ', 'Total Emission of BC               ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_emission'  )
          mixf(region)%f2d(emiso2      )%mf  = metafields( -1     , 'emiso2      ', 'Total Emission of SO2              ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_emission'  )
          mixf(region)%f2d(emiso4      )%mf  = metafields( -1     , 'emiso4      ', 'Total Direct Emission of SO4       ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_emission'  )
          mixf(region)%f2d(emidust     )%mf  = metafields( -1     , 'emidust     ', 'Total Emission of Dust             ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_emission'  )
          mixf(region)%f2d(emidms      )%mf  = metafields( -1     , 'emidms      ', 'Total Emission of DMS              ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_emission'  )
          mixf(region)%f2d(emiss       )%mf  = metafields( -1     , 'emiss       ', 'Total Emission of SeaSalt          ' , 'kg m-2 s-1', 'up',   'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_emission'  )
          mixf(region)%f2d(dryso2      )%mf  = metafields( -1     , 'dryso2      ', 'Dry Deposition of SO2              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_dry_deposition'  )
          mixf(region)%f2d(dryoa       )%mf  = metafields( -1     , 'dryoa       ', 'Dry Deposition of POM              ' , 'kg m-2 s-1', 'down', &
                           'tendency_of_atmosphere_mass_content_of_primary_particulate_organic_matter_dry_aerosol_due_to_dry_deposition'  )
          mixf(region)%f2d(drybc       )%mf  = metafields( -1     , 'drybc       ', 'Dry Deposition of BC               ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_dry_deposition'  )
          mixf(region)%f2d(dryso4      )%mf  = metafields( -1     , 'dryso4      ', 'Dry Deposition of SO4              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_dry_deposition'  )
          mixf(region)%f2d(drydust     )%mf  = metafields( -1     , 'drydust     ', 'Dry Deposition of Dust             ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_dry_deposition'  )
          mixf(region)%f2d(drydms      )%mf  = metafields( -1     , 'drydms      ', 'Dry Deposition of DMS              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_dry_deposition'  )
          mixf(region)%f2d(dryss       )%mf  = metafields( -1     , 'dryss       ', 'Dry Deposition of SeaSalt          ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_dry_deposition'  )
          mixf(region)%f2d(wetoa       )%mf  = metafields( -1     , 'wetoa       ', 'Wet Deposition of POM              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_particulate_organic_matter_dry_aerosol_due_to_wet_deposition'  )
          mixf(region)%f2d(wetbc       )%mf  = metafields( -1     , 'wetbc       ', 'Wet Deposition of BC               ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_black_carbon_dry_aerosol_due_to_wet_deposition'  )
          mixf(region)%f2d(wetso2      )%mf  = metafields( -1     , 'wetso2      ', 'Wet Deposition of SO2              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfur_dioxide_due_to_wet_deposition'  )
          mixf(region)%f2d(wetso4      )%mf  = metafields( -1     , 'wetso4      ', 'Wet Deposition of SO4              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_sulfate_dry_aerosol_due_to_wet_deposition'  )
          mixf(region)%f2d(wetdust     )%mf  = metafields( -1     , 'wetdust     ', 'Wet Deposition of Dust             ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dust_dry_aerosol_due_to_wet_deposition'  )
          mixf(region)%f2d(wetdms      )%mf  = metafields( -1     , 'wetdms      ', 'Wet Deposition of DMS              ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_dimethyl_sulfide_due_to_wet_deposition'  )
          mixf(region)%f2d(wetss       )%mf  = metafields( -1     , 'wetss       ', 'Wet Deposition of Seasalt          ' , 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_seasalt_dry_aerosol_due_to_wet_deposition'  )
          mixf(region)%f2d(od550aer    )%mf  = metafields( -1     , 'od550aer    ', 'AOD@550nm                          ' , '1         ', '',     'atmosphere_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(od550so4    )%mf  = metafields( -1     , 'od550so4    ', 'SO4 AOD@550nm                      ' , '1         ', '',     'atmosphere_optical_thickness_due_to_sulfate_ambient_aerosol'  )
          mixf(region)%f2d(od550bc     )%mf  = metafields( -1     , 'od550bc     ', 'Black Carbon AOD@550nm             ' , '1         ', '',     'atmosphere_optical_thickness_due_to_black_carbon_ambient_aerosol'  )
          mixf(region)%f2d(od550oa     )%mf  = metafields( -1     , 'od550oa     ', 'POM AOD@550nm                      ' , '1         ', '',     'atmosphere_optical_thickness_due_to_particulate_organic_matter_ambient_aerosol'  )
          mixf(region)%f2d(od550ss     )%mf  = metafields( -1     , 'od550ss     ', 'SeaSalt AOD@550nm                  ' , '1         ', '',     'atmosphere_optical_thickness_due_to_seasalt_ambient_aerosol'  )
          mixf(region)%f2d(od550dust   )%mf  = metafields( -1     , 'od550dust   ', 'Dust AOD@550nm                     ' , '1         ', '',     'atmosphere_optical_thickness_due_to_dust_ambient_aerosol'  )
          mixf(region)%f2d(od550no3    )%mf  = metafields( -1     , 'od550no3    ', 'Nitrate AOD@550nm                  ' , '1         ', '',     'atmosphere_optical_thickness_due_to_nitrate_ambient_aerosol'  )
          mixf(region)%f2d(od550aerh2o )%mf  = metafields( -1     , 'od550aerh2o ', 'Aerosol Water AOD@550nm            ' , '1         ', '',     'atmosphere_optical_thickness_due_to_water_in_ambient_aerosol'  )
          mixf(region)%f2d(od550lt1aer )%mf  = metafields( -1     , 'od550lt1aer ', 'Fine Mode AOD@550nm                ' , '1         ', '',     'atmosphere_optical_thickness_due_to_pm1_ambient_aerosol'  )
          mixf(region)%f2d(od550lt1dust)%mf  = metafields( -1     , 'od550lt1dust', 'Fine Mode Dust AOD@550nm           ' , '1         ', '',     'atmosphere_optical_thickness_due_to_dust_pm1_ambient_aerosol'  )
          mixf(region)%f2d(od550lt1ss)%mf    = metafields( -1     , 'od550lt1ss  ', 'Fine Mode SeaSalt AOD@550nm        ' , '1         ', '',     'atmosphere_optical_thickness_due_to_seasalt_pm1_ambient_aerosol'  )
          mixf(region)%f2d(abs550aer   )%mf  = metafields( -1     , 'abs550aer   ', 'Absorption AOD@550nm               ' , '1         ', '',     'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(asyaer      )%mf  = metafields( -1     , 'asyaer      ', 'Asymmetry Parameter                ' , '1         ', '',     'atmosphere_asymmetry_parameter_ambient_aerosol'  )
          mixf(region)%f2d(ec550aer    )%mf  = metafields( -1     , 'ec550aer    ', 'Surface Ambient Aerosol Extinction@550nm' , 'm-1  ', '',     'surface_extinction_due_to_ambient_aerosol'  )
          mixf(region)%f2d(ec550dryaer )%mf  = metafields( -1     , 'ec550dryaer ', 'Surface Dry Aerosol Extinction@550 nm','m-1       ', '',     'surface_extinction_due_to_dry_aerosol'  )
          mixf(region)%f2d(abs550dryaer)%mf  = metafields( -1     , 'abs550dryaer', 'Surface Dry Aerosol Absorption@550 nm','m-1       ', '',     'surface_absorption_due_to_dry_aerosol'  )
          mixf(region)%f2d(asydryaer   )%mf  = metafields( -1     , 'asydryaer   ', 'Surface Dry Aersol Asymmetry Parameter' , '1      ', '',     'surface_asymmetry_parameter_dry_aerosol'  )
          mixf(region)%f2d(ec550drylt1aer )%mf  = metafields( -1  , 'ec550drylt1aer ', 'Surface Fine Mode Dry Aerosol Extinction@550 nm' , 'm-1', '', 'surface_extinction_due_to_pm1_dry_aerosol'  )
          mixf(region)%f2d(abs550drylt1aer)%mf  = metafields( -1  , 'abs550drylt1aer', 'Surface Fine Mode Dry Aerosol Absorption@550 nm' , 'm-1', '', 'surface_absorption_due_to_pm1_dry_aerosol'  )
          mixf(region)%f2d(od440aer    )%mf  = metafields( -1     , 'od440aer    ', 'AOD@440nm                          ' , '1         ', '',     'atmosphere_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(abs440aer   )%mf  = metafields( -1     , 'abs440aer   ', 'Absorption AOD@440nm               ' , '1         ', '',     'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(ec440dryaer )%mf  = metafields( -1     , 'ec440dryaer ', 'Surface Dry Aerosol Extinction@440nm' , 'm-1      ', '',     'surface_extinction_due_to_dry_aerosol'  )
          mixf(region)%f2d(abs440dryaer)%mf  = metafields( -1     , 'abs440dryaer', 'Surface Dry Aerosol Absorption@440nm' , 'm-1      ', '',     'surface_absorption_due_to_dry_aerosol'  )
          mixf(region)%f2d(od870aer    )%mf  = metafields( -1     , 'od870aer    ', 'AOD@870nm                          ' , '1         ', '',     'atmosphere_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(abs870aer   )%mf  = metafields( -1     , 'abs870aer   ', 'Absorption AOD@870nm               ' , '1         ', '',     'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(ec870dryaer )%mf  = metafields( -1     , 'ec870dryaer ', 'Surface Dry Aerosol Extinction@870nm' , 'm-1      ', '',     'surface_extinction_due_to_dry_aerosol'  )
          mixf(region)%f2d(abs870dryaer)%mf  = metafields( -1     , 'abs870dryaer', 'Surface Dry Aerosol Absorption@870nm' , 'm-1      ', '',     'surface_absorption_due_to_dry_aerosol'  )
          mixf(region)%f2d(od350aer    )%mf  = metafields( -1     , 'od350aer    ', 'AOD@350nm                          ' , '1         ', '',     'atmosphere_optical_thickness_due_to_ambient_aerosol'  )
          mixf(region)%f2d(abs350aer   )%mf  = metafields( -1     , 'abs350aer   ', 'Absorption AOD@350nm               ' , '1         ', '',     'atmosphere_absorption_optical_thickness_due_to_ambient_aerosol'  )


!!$        mixf(region)%f2d(52)%mf  = metafields( 'od550soa    ', 'AOD@550nm SOA                      ' , '1         ' )
!!$        mixf(region)%f2d(53)%mf  = metafields( 'od550bb     ', 'AOD@550nm BB                       ' , '1         ' )
!!$        mixf(region)%f2d(54)%mf  = metafields( 'gf90aer     ', 'GF @ 90 % RH                       ' , '1         ' )

          if (stat_output) then
             ! 1d fields (stations)
             statf(region)%f1d( 1)%mf  = metafields2( iso4nus, 'mmr1Dtr01',      'mmr of tracer 01' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 2)%mf  = metafields2( iso4ais, 'mmr1Dtr02',      'mmr of tracer 02' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 3)%mf  = metafields2( iso4acs, 'mmr1Dtr03',      'mmr of tracer 03' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 4)%mf  = metafields2( iso4cos, 'mmr1Dtr04',      'mmr of tracer 04' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 5)%mf  = metafields2( ibcais , 'mmr1Dtr05',      'mmr of tracer 05' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 6)%mf  = metafields2( ibcacs , 'mmr1Dtr06',      'mmr of tracer 06' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 7)%mf  = metafields2( ibccos , 'mmr1Dtr07',      'mmr of tracer 07' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 8)%mf  = metafields2( ibcaii , 'mmr1Dtr08',      'mmr of tracer 08' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d( 9)%mf  = metafields2( ipomais, 'mmr1Dtr09',      'mmr of tracer 09' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(10)%mf  = metafields2( ipomacs, 'mmr1Dtr10',      'mmr of tracer 10' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(11)%mf  = metafields2( ipomcos, 'mmr1Dtr11',      'mmr of tracer 11' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(12)%mf  = metafields2( ipomaii, 'mmr1Dtr12',      'mmr of tracer 12' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(13)%mf  = metafields2( issacs , 'mmr1Dtr13',      'mmr of tracer 13' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(14)%mf  = metafields2( isscos , 'mmr1Dtr14',      'mmr of tracer 14' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(15)%mf  = metafields2( iduacs , 'mmr1Dtr15',      'mmr of tracer 15' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(16)%mf  = metafields2( iducos , 'mmr1Dtr16',      'mmr of tracer 16' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(17)%mf  = metafields2( iduaci , 'mmr1Dtr17',      'mmr of tracer 17' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(18)%mf  = metafields2( iducoi , 'mmr1Dtr18',      'mmr of tracer 18' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f1d(19)%mf  = metafields2( inus_n , 'conccn1Dmode01', 'number concentration of mode 01' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(20)%mf  = metafields2( iais_n , 'conccn1Dmode02', 'number concentration of mode 02' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(21)%mf  = metafields2( iacs_n , 'conccn1Dmode03', 'number concentration of mode 03' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(22)%mf  = metafields2( icos_n , 'conccn1Dmode04', 'number concentration of mode 04' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(23)%mf  = metafields2( iaii_n , 'conccn1Dmode05', 'number concentration of mode 05' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(24)%mf  = metafields2( iaci_n , 'conccn1Dmode06', 'number concentration of mode 06' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(25)%mf  = metafields2( icoi_n , 'conccn1Dmode07', 'number concentration of mode 07' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f1d(26)%mf  = metafields2(  1     , 'mmr1Daerh2o01 ', 'mmr of aerosol water of mode 01' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f1d(27)%mf  = metafields2(  2     , 'mmr1Daerh2o02 ', 'mmr of aerosol water of mode 02' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f1d(28)%mf  = metafields2(  3     , 'mmr1Daerh2o03 ', 'mmr of aerosol water of mode 03' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f1d(29)%mf  = metafields2(  4     , 'mmr1Daerh2o04 ', 'mmr of aerosol water of mode 04' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f1d(30)%mf  = metafields2( ino3_a , 'mmr1Dno3_a    ', 'mmr of nitrate aerosol' ,  'kg kg-1', 'mass_fraction_of_nitrate_dry_aerosol_in_air'  )
             statf(region)%f1d(31)%mf  = metafields2( inh4   , 'mmr1Dnh4      ', 'mmr of ammonium' ,  'kg kg-1', 'mass_fraction_of_ammonium_dry_aerosol_in_air'  )
             statf(region)%f1d(32)%mf  = metafields2( imsa   , 'mmr1Dmsa      ', 'mmr of methane sulfonic acid' ,  'kg kg-1', 'mass_fraction_of_methane_sulfonic_acid_in_air'  )
             statf(region)%f1d(33)%mf  = metafields2( -1     , 'ec550dry1Daer ', 'dry aerosol extinction at 550 nm' ,      'm-1'    , 'atmosphere_extinction_due_to_dry_aerosol'  )
             statf(region)%f1d(34)%mf  = metafields2( -1     , 'abs550dry1Daer', 'dry aerosol absorption at 550 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(35)%mf  = metafields2( -1     , 'asydry1Daer', 'dry aerosol asymmetry parameter' ,      '1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(36)%mf  = metafields2( -1     , 'ec550dry1Dlt1aer', 'fine mode dry aerosol extinction at 550 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(37)%mf  = metafields2( -1     , 'abs550dry1Dlt1aer', 'fine mode dry aerosol absorption at 550 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(38)%mf  = metafields2( -1     , 'ec440dry1Daer', 'dry aerosol extinction at 440 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(39)%mf  = metafields2( -1     , 'abs440dry1Daer', 'dry aerosol absorption at 440 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(40)%mf  = metafields2( -1     , 'ec870dry1Daer', 'dry aerosol extinction at 870 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(41)%mf  = metafields2( -1     , 'abs870dry1Daer', 'dry aerosol absorption at 870 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f1d(42)%mf  = metafields2( -1     , 'temp1D          ', 'Temperature                     ' ,      'K'      , 'air_temperature'  )
             statf(region)%f1d(43)%mf  = metafields2( -1     , 'hus1D           ', 'Specific Humidity               ' ,      '1'      , 'specific_humidity'  )
             statf(region)%f1d(44)%mf  = metafields2( -1     , 'pressure1D      ', 'Pressure                        ' ,      'Pa'     , 'air_pressure'  )

             ! 0d fields (stations)
             statf(region)%f0d( 1)%mf  = metafields2( iso4nus, 'mmrtr01',      'mmr of tracer 01' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 2)%mf  = metafields2( iso4ais, 'mmrtr02',      'mmr of tracer 02' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 3)%mf  = metafields2( iso4acs, 'mmrtr03',      'mmr of tracer 03' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 4)%mf  = metafields2( iso4cos, 'mmrtr04',      'mmr of tracer 04' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 5)%mf  = metafields2( ibcais , 'mmrtr05',      'mmr of tracer 05' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 6)%mf  = metafields2( ibcacs , 'mmrtr06',      'mmr of tracer 06' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 7)%mf  = metafields2( ibccos , 'mmrtr07',      'mmr of tracer 07' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 8)%mf  = metafields2( ibcaii , 'mmrtr08',      'mmr of tracer 08' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d( 9)%mf  = metafields2( ipomais, 'mmrtr09',      'mmr of tracer 09' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(10)%mf  = metafields2( ipomacs, 'mmrtr10',      'mmr of tracer 10' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(11)%mf  = metafields2( ipomcos, 'mmrtr11',      'mmr of tracer 11' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(12)%mf  = metafields2( ipomaii, 'mmrtr12',      'mmr of tracer 12' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(13)%mf  = metafields2( issacs , 'mmrtr13',      'mmr of tracer 13' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(14)%mf  = metafields2( isscos , 'mmrtr14',      'mmr of tracer 14' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(15)%mf  = metafields2( iduacs , 'mmrtr15',      'mmr of tracer 15' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(16)%mf  = metafields2( iducos , 'mmrtr16',      'mmr of tracer 16' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(17)%mf  = metafields2( iduaci , 'mmrtr17',      'mmr of tracer 17' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(18)%mf  = metafields2( iducoi , 'mmrtr18',      'mmr of tracer 18' ,                 'kg kg-1', 'mass_fraction_of_tracer_dry_aerosol_in_air'  )
             statf(region)%f0d(19)%mf  = metafields2( inus_n , 'conccnmode01', 'number concentration of mode 01' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(20)%mf  = metafields2( iais_n , 'conccnmode02', 'number concentration of mode 02' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(21)%mf  = metafields2( iacs_n , 'conccnmode03', 'number concentration of mode 03' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(22)%mf  = metafields2( icos_n , 'conccnmode04', 'number concentration of mode 04' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(23)%mf  = metafields2( iaii_n , 'conccnmode05', 'number concentration of mode 05' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(24)%mf  = metafields2( iaci_n , 'conccnmode06', 'number concentration of mode 06' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(25)%mf  = metafields2( icoi_n , 'conccnmode07', 'number concentration of mode 07' ,  'm-3'    , 'number_concentration_of_ambient_aerosol_in_air'  )
             statf(region)%f0d(26)%mf  = metafields2(  1     , 'mmraerh2o01 ', 'mmr of aerosol water of mode 01' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f0d(27)%mf  = metafields2(  2     , 'mmraerh2o02 ', 'mmr of aerosol water of mode 02' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f0d(28)%mf  = metafields2(  3     , 'mmraerh2o03 ', 'mmr of aerosol water of mode 03' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f0d(29)%mf  = metafields2(  4     , 'mmraerh2o04 ', 'mmr of aerosol water of mode 04' ,  'kg kg-1', 'mass_fraction_of_water_in_ambient_aerosol_in_air'  )
             statf(region)%f0d(30)%mf  = metafields2( ino3_a , 'mmrno3_a    ', 'mmr of nitrate aerosol' ,  'kg kg-1', 'mass_fraction_of_nitrate_dry_aerosol_in_air'  )
             statf(region)%f0d(31)%mf  = metafields2( inh4   , 'mmrnh4      ', 'mmr of ammonium' ,  'kg kg-1', 'mass_fraction_of_ammonium_dry_aerosol_in_air'  )
             statf(region)%f0d(32)%mf  = metafields2( imsa   , 'mmrmsa      ', 'mmr of methane sulfonic acid' ,  'kg kg-1', 'mass_fraction_of_methane_sulfonic_acid_in_air'  )
             statf(region)%f0d(33)%mf  = metafields2( -1     , 'ec550dryaer ', 'dry aerosol extinction at 550 nm' ,      'm-1'    , 'atmosphere_extinction_due_to_dry_aerosol'  )
             statf(region)%f0d(34)%mf  = metafields2( -1     , 'abs550dryaer', 'dry aerosol absorption at 550 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(35)%mf  = metafields2( -1     , 'asydryaer', 'dry aerosol asymmetry parameter' ,      '1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(36)%mf  = metafields2( -1     , 'ec550drylt1aer', 'fine mode dry aerosol extinction at 550 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(37)%mf  = metafields2( -1     , 'abs550drylt1aer', 'fine mode dry aerosol absorption at 550 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(38)%mf  = metafields2( -1     , 'ec440dryaer', 'dry aerosol extinction at 440 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(39)%mf  = metafields2( -1     , 'abs440dryaer', 'dry aerosol absorption at 440 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(40)%mf  = metafields2( -1     , 'ec870dryaer', 'dry aerosol extinction at 870 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(41)%mf  = metafields2( -1     , 'abs870dryaer', 'dry aerosol absorption at 870 nm' ,      'm-1'    , 'atmosphere_absorption_due_to_dry_aerosol'  )
             statf(region)%f0d(42)%mf  = metafields2( -1     , 'temp        ', 'Temperature                     ' ,      'K'      , 'air_temperature'  )
             statf(region)%f0d(43)%mf  = metafields2( -1     , 'hus         ', 'Specific Humidity               ' ,      '1'      , 'specific_humidity'  )
             statf(region)%f0d(44)%mf  = metafields2( -1     , 'pressure    ', 'Pressure                        ' ,      'Pa'     , 'air_pressure'  ) 
             statf(region)%f0d(45)%mf  = metafields2( -1     , 'ps          ', 'Surface Air Pressure            ' ,      'Pa '    , 'surface_air_pressure'  )
             statf(region)%f0d(46)%mf  = metafields2( -1     , 'precip      ', 'Precipitation                   ' ,      'kg m-2 s-1', 'precipitation_flux'  )
          endif

          ! set global dimensions of fields for netcdf definitions
          mixf (region)%nlon      = imr
          mixf (region)%nlat      = jmr
          mixf (region)%nlev      = lmr
          if (stat_output) then 
             statf(region)%nstat     = nstations
             statf(region)%nlev      = lmr
          endif

          ! intermediate structures for budgets
          allocate(drydepos(region)%f2dslast(i1:i2,j1:j2,7))
          allocate(wetdepos(region)%f2dslast(i1:i2,j1:j2,7))
          allocate(emission(region)%f2dslast(i1:i2,j1:j2,7))

          ! these here are the initial budgets (monthly): 0.0
          drydepos(region)%f2dslast = 0.0
          wetdepos(region)%f2dslast = 0.0
          emission(region)%f2dslast = 0.0

       endif ! newsrun

       ! reset counters and accumulators
       mixf (region)%acct = 0
       do i = 1, ntracer_3d
          mixf(region)%f3d(i)%field = 0.0
       end do
       do i = 1, ntracer_2d
          mixf(region)%f2d(i)%field = 0.0
       end do
       if (stat_output) then
          ! stations
          statf(region)%acct = 0
          do i = 1, ntracer_1dstat
             statf(region)%f1d(i)%field = 0.0
          end do
          do i = 1, ntracer_0dstat
             statf(region)%f0d(i)%field = 0.0
          end do
       endif

       ! ----------------
       ! open NetCDF file  (mixf)
       ! ----------------

       call tau2date(itau,idater)

       ! time (in days since 2001-01-01 00:00) 
       call date2tau( time_reftime6, itauref )

       select case (trim(aerocom_freq)) 
         case ('hourly')
           write(idates,  '(i4,3i2.2)') idater(1), idater(2), idater(3), idater(4)+1
           write(lidates, '(i4,"-",i2.2,"-",i2.2," ",i2.2,":30")') idater(1), idater(2), idater(3), idater(4)
           idater(5:6) = (/30,0/) 
         case ('daily')
           write(idates,  '(i4,2i2.2)') idater(1), idater(2), idater(3)
           write(lidates, '(i4,"-",i2.2,"-",i2.2," 12:00")') idater(1), idater(2), idater(3)
           ! set noon (recommendations)
           idater(4:6) = (/12,0,0/)
         case ('monthly')
           ! for monthly files, set time to middle of the month
           write(idates,  '(i4,i2.2)') idater(1), idater(2)
           mlength=mlen(idater(2))
           if ( mod(mlength,2) == 0 ) then
             idater(3:6) = (/mlength/2 + 1,0,0,0/)
             write(lidates, '(i4,"-",i2.2,"-",i2.2," 00:00")') idater(1), idater(2), idater(3)
           else
             idater(3:6) = (/(mlength+1)/2,12,0,0/)
             write(lidates, '(i4,"-",i2.2,"-",i2.2," 12:00")') idater(1), idater(2), idater(3)
           endif
         case default
           write (gol,'("Unknown AeroCom output frequency;")'); call goErr
       end select

       call date2tau( idater, itaucur )
       reftime = (itaucur - itauref) / 86400.
       
       ! Changed file name convention to AeroCom Control 2015
       mixf(region)%fname = trim(outdir)//'/'//&
            trim(f_dataid)//'_'//&
            trim(f_model) //'_'//&
            trim(aerocom_exper)//'_'//&
            trim(f_dimmix)//'_'//&
            trim(idates)  //'_'//&
            trim(aerocom_freq)  //'.nc'

#ifdef MPI     
       ! overwrite existing files (clobber), provide MPI stuff:
       call MDF_Create( mixf(region)%fname, MDF_NETCDF4, MDF_REPLACE, mixf(region)%funit, 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( mixf(region)%fname, MDF_NETCDF4, MDF_REPLACE, mixf(region)%funit, status )
       IF_NOTOK_RETURN(status=1)
#endif

       if(okdebug) then
          write(gol,*)  'output_aerocom_init: File ', trim(mixf(region)%fname), ' opened ' ; call goPr
       endif

       ! global attributes
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'title',               'Model output for AeroCom', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'source',              'TM5-mp', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'institution',         'Royal Netherlands Meteorological Institute (KNMI), De Bilt, The Netherlands)', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'contact'     ,        'Twan van Noije; noije@knmi.nl', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'project_id',          'AeroCom Phase 3', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'conventions',         'CF-1.0 - HTAP', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'date',                 lidates, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'time_units',           time_units, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'references',          'http://tm5.sourceforge.net/', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! define dimensions
       call MDF_Def_Dim( mixf(region)%funit, 'longitude',     imr, mixf(region)%nlon, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Def_Dim( mixf(region)%funit, 'latitude',      jmr, mixf(region)%nlat, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Def_Dim( mixf(region)%funit, 'alevel',        lmr, mixf(region)%nlev, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! define variables
       do i = 1, ntracer_3d 
          call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f3d(i)%mf%vname), MDF_FLOAT, &
               (/mixf(region)%nlon, mixf(region)%nlat, mixf(region)%nlev/), varid, status )
          IF_NOTOK_MDF(fid=mixf(region)%funit) 
          call MDF_Var_Par_Access( mixf(region)%funit, varid, access_mode, status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'long_name',     trim(mixf(region)%f3d(i)%mf%lname), status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'standard_name', trim(mixf(region)%f3d(i)%mf%standard_name), status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'units',         trim(mixf(region)%f3d(i)%mf%unit), status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'time',          reftime, status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'cell_methods',  'time: mean', status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          if( mixf(region)%f3d(i)%mf%positive /= '' ) then 
             call MDF_Put_Att( mixf(region)%funit, varid, 'positive',  trim(mixf(region)%f3d(i)%mf%positive), status )
             IF_NOTOK_MDF(fid=mixf(region)%funit)
          end if
          mixf(region)%f3d(i)%varid = varid
       end do
       do i = 1, ntracer_2d 
          call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f2d(i)%mf%vname), MDF_FLOAT, &
               (/mixf(region)%nlon, mixf(region)%nlat/), varid, status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Var_Par_Access( mixf(region)%funit, varid, access_mode, status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'long_name',     trim(mixf(region)%f2d(i)%mf%lname), status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'standard_name', trim(mixf(region)%f2d(i)%mf%standard_name), status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'units',         trim(mixf(region)%f2d(i)%mf%unit), status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'time',          reftime, status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          call MDF_Put_Att( mixf(region)%funit, varid, 'cell_methods',  'time: mean', status )
          IF_NOTOK_MDF(fid=mixf(region)%funit)
          if( mixf(region)%f2d(i)%mf%positive /= '' ) then 
             call MDF_Put_Att( mixf(region)%funit, varid, 'positive',  trim(mixf(region)%f2d(i)%mf%positive), status )
             IF_NOTOK_MDF(fid=mixf(region)%funit)
          end if
          mixf(region)%f2d(i)%varid = varid
       end do

       call MDF_EndDef( mixf(region)%funit, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       if (stat_output) then
         ! ------------
         ! station NetCDF file 
         ! ------------
         statf(region)%fname = trim(outdir)//'/'//&
               trim(f_dataid) //'_'//&
               trim(f_model)  //'_'//&
               trim(aerocom_exper) //'_'//&
               trim(f_dimstat)//'_'//&
               trim(idates)   //'_'//&
               trim(aerocom_freq)   //'.nc'

#ifdef MPI
          ! overwrite existing files (clobber), provide MPI stuff:
          call MDF_Create( statf(region)%fname, MDF_NETCDF4, MDF_REPLACE, statf(region)%funit, 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( statf(region)%fname, MDF_NETCDF4, MDF_REPLACE, statf(region)%funit, status )
          IF_NOTOK_RETURN(status=1)
#endif

          if(okdebug) then
             write(gol,*) ' output_aerocom_init: File ', trim(statf(region)%fname), ' opened '; call goPr
          endif

          ! global attributes
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'title',       'Model output for AeroCom', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'source',      'TM5-mp', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'institution', 'Royal Netherlands Meteorological Institute (KNMI), De Bilt, The Netherlands)', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'contact',     'Twan van Noije; noije@knmi.nl', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'project_id',  'AeroCom Phase 3', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'conventions', 'CF-1.0 - HTAP', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'date',         lidates, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'time_units',   time_units, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, MDF_GLOBAL, 'references',  'http://tm5.sourceforge.net/', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)

          ! define dimensions
          call MDF_Def_Dim( statf(region)%funit, 'station',         nstations,    statf(region)%nstat, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Def_Dim( statf(region)%funit, 'stationname_len', statnamelen,  snamelendim, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Def_Dim( statf(region)%funit, 'alevel',          lmr,          statf(region)%nlev, status )
          IF_NOTOK_MDF(fid=statf(region)%funit) 

          ! define variables
          do i = 1, ntracer_1dstat
             call MDF_Def_Var( statf(region)%funit, trim(statf(region)%f1d(i)%mf%vname), MDF_FLOAT, &
                  (/statf(region)%nstat,statf(region)%nlev/), varid, status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'long_name',     trim(statf(region)%f1d(i)%mf%lname), status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', trim(statf(region)%f1d(i)%mf%standard_name), status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'units',         trim(statf(region)%f1d(i)%mf%unit), status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'time',          reftime, status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'cell_methods',  'time: mean', status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             statf(region)%f1d(i)%varid = varid
          end do

          do i = 1, ntracer_0dstat 
             call MDF_Def_Var( statf(region)%funit, trim(statf(region)%f0d(i)%mf%vname), MDF_FLOAT, &
                  (/statf(region)%nstat/), varid, status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'long_name',     trim(statf(region)%f0d(i)%mf%lname), status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'standard_name', trim(statf(region)%f0d(i)%mf%standard_name), status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'units',         trim(statf(region)%f0d(i)%mf%unit), status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'time',          reftime, status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             call MDF_Put_Att( statf(region)%funit, varid, 'cell_methods',  'time: mean', status )
             IF_NOTOK_MDF(fid=statf(region)%funit)
             statf(region)%f0d(i)%varid = varid
          end do


          ! station information (name, nb, lat, lon, alt)
          call MDF_Def_Var( statf(region)%funit, 'stationname', MDF_CHAR, (/snamelendim,statf(region)%nstat/), varid, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
          IF_NOTOK_MDF(fid=statf(region)%funit) 
          call MDF_Put_Att( statf(region)%funit, varid, 'long_name',      'station name', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'standard_name',  'station_name', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'units',          '1', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          statvarid(1) = varid 

          call MDF_Def_Var( statf(region)%funit, 'stationnb', MDF_INT, (/statf(region)%nstat/), varid, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'long_name',      'station number', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'standard_name',  'station_nb', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'units',          '1', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          statvarid(2) = varid 

          call MDF_Def_Var( statf(region)%funit, 'stationlat', MDF_FLOAT, (/statf(region)%nstat/), varid, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'long_name',      'location of station: station latitude', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'standard_name',  'stationlat', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'units',          '1', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          statvarid(3) = varid 

          call MDF_Def_Var( statf(region)%funit, 'stationlon', MDF_FLOAT, (/statf(region)%nstat/), varid, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'long_name',      'location of station: station longitude', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'standard_name',  'stationlon', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'units',          '1', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          statvarid(4) = varid 

          call MDF_Def_Var( statf(region)%funit, 'stationalt', MDF_FLOAT, (/statf(region)%nstat/), varid, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Var_Par_Access( statf(region)%funit, varid, access_mode_sta, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'long_name',      'location of station: station altitude', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'standard_name',  'stationalt', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Att( statf(region)%funit, varid, 'units',          '1', status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          statvarid(5) = varid 

          call MDF_EndDef( statf(region)%funit, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)

          ! put station information (note: behavior here is like if all procs write same data)
          call MDF_Put_Var( statf(region)%funit, statvarid(1), tmpstatnames, status, start=(/ 1, 1/), count=(/ statnamelen,nstations /) )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Var( statf(region)%funit, statvarid(2), statlist(:)%statnb, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Var( statf(region)%funit, statvarid(3), statlist(:)%statlat, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Var( statf(region)%funit, statvarid(4), statlist(:)%statlon, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
          call MDF_Put_Var( statf(region)%funit, statvarid(5), statlist(:)%statalt, status )
          IF_NOTOK_MDF(fid=statf(region)%funit)
       endif

    end do regionloop ! nregions

    call goLabel() ; status = 0

  end  subroutine output_aerocom_init
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_AEROCOM_DONE
  !
  ! !DESCRIPTION: Free parameters. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_aerocom_done(stat_output, status, iregion)
    !
    ! !USES:
    !
    ! !INPUT PARAMETERS:
    !
    logical, intent(in)           :: stat_output    ! true if stations output is requested
    integer, intent(in), optional :: iregion
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)           :: status
    !
    ! !REVISION HISTORY: 
    !   29 Nov 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/output_aerocom_done'
    
    integer :: i, region
    
    ! --- begin -------------------------------------

    call goLabel(rname)

    deallocate( wdep_out )

    regionloop: do region = 1, nregions

       ! if region given, cycle if other region!
       if (present(iregion)) then
          if(iregion /= region) cycle regionloop
       endif

       do i=1,ntracer_3d
          deallocate( mixf(region)%f3d(i)%field )
       end do
       do i=1,ntracer_2d
          deallocate( mixf(region)%f2d(i)%field )
       end do
       deallocate( mixf (region)%f3d )
       deallocate( mixf (region)%f2d )
       deallocate( drydepos(region)%f2dslast )
       deallocate( wetdepos(region)%f2dslast )
       deallocate( emission(region)%f2dslast )

       if (stat_output) then
          ! stations
          do i=1,ntracer_1dstat
             deallocate( statf(region)%f1d(i)%field )
          end do
          do i=1,ntracer_0dstat 
             deallocate( statf(region)%f0d(i)%field )
          end do
          deallocate( statf(region)%f1d )
          deallocate( statf(region)%f0d )
          deallocate( sintp(region)%ixl, sintp(region)%wixl, sintp(region)%ixr, sintp(region)%wixr, &
                      sintp(region)%jyl, sintp(region)%wjyl, sintp(region)%jyr, sintp(region)%wjyr, &
                      sintp(region)%lzl, sintp(region)%wlzl, sintp(region)%lzr, sintp(region)%wlzr, &
                      sintp(region)%tght )
       endif

    end do regionloop

    if (stat_output) then
       deallocate( statlist )
       deallocate( tmpstatnames )
       deallocate( stat_active )
    endif

    call goLabel() ; status = 0

  end subroutine output_aerocom_done
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  output_aerocom_write
  !
  ! !DESCRIPTION: This routines builds the average by dividing accumulated
  !               data by stack counter. The results are written to file. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_aerocom_write(region, stat_output, status)
    !
    ! !USES:
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    logical, intent(in)  :: stat_output    ! true if stations output is requested
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   29 Nov 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer              :: i, imr, jmr, lmr
    integer              :: i1, i2, j1, j2, ilev
    integer              :: istat

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

    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

    !  this here is already accumulated over the time period (day)
    call collect_budgets( region, status )
    IF_NOTOK_RETURN(status=1)

    ! ---------------------
    ! divide by accumulator
    ! ---------------------
    do i = 1, ntracer_3d
       mixf(region)%f3d(i)%field = mixf(region)%f3d(i)%field / real( mixf(region)%acct )
    end do
    do i = 1, ntracer_2d
       mixf(region)%f2d(i)%field = mixf(region)%f2d(i)%field / real( mixf(region)%acct )
    end do
    if (stat_output) then
       ! stations
       do istat = 1, nstations
          if (stat_active(istat)) then
             do i = 1, ntracer_1dstat
                statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) / real( statf(region)%acct )
             end do
             do i = 1, ntracer_0dstat
                statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) / real( statf(region)%acct )
             end do
          endif
       end do
    endif

    select case (trim(aerocom_freq))
      case ('hourly')
        write(gol,'("---> AEROCOM diagnostics: write file for previous hour")'); call goPr
      case ('daily')
        write(gol,'("---> AEROCOM diagnostics: write file for previous day")'); call goPr
      case ('monthly') 
        write(gol,'("---> AEROCOM diagnostics: write file for previous month")'); call goPr
    end select

    ! -------------
    ! write to file 
    ! -------------
    do i=1,ntracer_2d
       call MDF_Put_Var( mixf(region)%funit, mixf(region)%f2d(i)%varid, mixf(region)%f2d(i)%field, status, start=(/i1,j1/), count=(/imr,jmr/) )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
    end do
    do i=1,ntracer_3d
       call MDF_Put_Var( mixf(region)%funit, mixf(region)%f3d(i)%varid, mixf(region)%f3d(i)%field, status, start=(/i1,j1,1/), count=(/imr,jmr,lmr/) )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
    end do
    call MDF_Close( mixf(region)%funit, status )
    IF_NOTOK_MDF(fid=mixf(region)%funit)

    if (stat_output) then
       ! stations
       do istat = 1, nstations
          if (stat_active(istat)) then
             do i=1,ntracer_1dstat
                call MDF_Put_Var( statf(region)%funit, statf(region)%f1d(i)%varid, statf(region)%f1d(i)%field(istat,:), &
                     status, start=(/istat,1/), count=(/1,lmr/) )
                IF_NOTOK_MDF(fid=statf(region)%funit)
             end do
             do i=1,ntracer_0dstat
                call MDF_Put_Var( statf(region)%funit, statf(region)%f0d(i)%varid, (/statf(region)%f0d(i)%field(istat)/), &
                     status, start=(/istat/), count=(/1/) )
                IF_NOTOK_MDF(fid=statf(region)%funit)
             end do
          endif
       end do
       call MDF_Close( statf(region)%funit, status )
       IF_NOTOK_MDF(fid=statf(region)%funit)
    endif

    call goLabel() ; status = 0

  end subroutine output_aerocom_write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_AEROCOM_STEP
  !
  ! !DESCRIPTION: This is the accumulation routine. It is called following 
  !               the user specification aerocom.dhour in rc-file. It 
  !               evaluates the various diagnostics and does summing. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_aerocom_step( region, dhour, stat_output, status )
    !
    ! !USES:
    !
    use GO         , only : TDate, NewDate, rTotal, operator(-)
    use Grid       , only : FPressure
    use binas,       only : rgas, rol
    use datetime,    only : tau2date
    use MeteoData,   only : sp_dat, lsp_dat, cp_dat, m_dat
    use MeteoData,   only : temper_dat, humid_dat, gph_dat
    use global_data, only : mass_dat, region_dat
    use tracer_data, only : chem_dat
    use dims,        only : itaur
    use chem_param,  only : iso4acs, iso4cos, iduacs, iso4ais, ibccos, ibcaii, xmair
    use chem_param,  only : iso4nus, isscos, ino3_a, issacs,  iducos, iduaci, nmod
    use chem_param,  only : iducoi, ibcacs, ipomcos, ipomaii, ibcais, ipomacs, ipomais
    use chem_param,  only : imsa, inh4
    use chem_param,  only : ntracet
    use mo_aero_m7,  only : nsol
    use optics,      only : optics_aop_get
    use m7_data,     only : h2o_mode, rw_mode, rwd_mode
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    integer, intent(in)  :: dhour   ! this is aerocom.dhour [hours]
    logical, intent(in)  :: stat_output    ! true if stations output is requested
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   29 Nov 2010 - Achim Strunk - Creation
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC


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

    ! MPI arrays to gather fields.
    real, dimension(:,:,:,:),              pointer     :: rm, rmc
    real, dimension(:,:,:),   allocatable, target      :: mg
    real, dimension(:),                    pointer     :: dxyp

    real, dimension(:,:,:),                allocatable :: pres3d

    integer  :: i, j, k, imr, jmr, lmr, lwl, dtime
    integer  :: i1, i2, j1, j2
    integer  :: ie, je ! indices for subdomain extended with halo cells
    integer, parameter :: l_halo=1
    logical  :: polar
    integer  :: istat, imode
    real     :: dens, load_tmp

    Real, Dimension(:,:,:), Allocatable :: aop_output_ext ! extinctions
    Real, Dimension(:,:),   Allocatable :: aop_output_a   ! single scattering albedo
    Real, Dimension(:,:),   Allocatable :: aop_output_g   ! assymetry factor
    Real, Dimension(:,:,:), Allocatable :: aop_output_add ! additional parameters

    real, dimension(:,:,:,:,:), allocatable :: opt_ext
    real, dimension(:,:,:,:),   allocatable :: opt_a
    real, dimension(:,:,:,:),   allocatable :: opt_g
    real, dimension(:,:,:,:,:), allocatable :: opt_add
    real, dimension(:,:,:),     allocatable :: volume

    real, dimension(:,:),       allocatable :: temp2d

    integer :: ncontr, lvec, lct, l, indoffset, nwl
    integer :: nadd, nadd_max, indoffset_stat, indabs
    integer :: iadd
    real    :: no3fac, wght, dwght, tabs

    integer,dimension(6) :: idate_f

    type(TDate)          ::  t, t0
    real                 ::  time 

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

    call goLabel(rname)

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

    nullify(rmc)
    rmc    => chem_dat(region)%rm
    dxyp   => region_dat(region)%dxyp

    rm => mass_dat(region)%rm 

    ! get helping pressure field in 3d
    ie=i2; je=j2
    if (stat_output) then
      if (stat_halo_glob.eq.1) then
        CALL UPDATE_HALO( dgrid(region), sp_dat(region)%data(:,:,1), sp_dat(region)%halo, status)
        IF_NOTOK_RETURN(status=1)
        ! in the interpolation routine,
        ! halo cells are only needed at the eastern or northern boundary of the subdomain
        ie=i2+1
        if (.not.polar) je=j2+1
      endif
    endif
    allocate( pres3d(i1:ie,j1:je,lmr) )
    ! fill mid level pressure
    call FPressure( levi, sp_dat(region)%data(i1:ie,j1:je,1), pres3d, status )
    IF_NOTOK_RETURN(status=1)

    ! dtime is seconds per step 
    ! dtime will be taken as weight for summing up 
    ! (makes it fit for varying step lengths) 
    dtime = dhour * 3600. 

    ! sum up for later averaging
    mixf (region)%acct = mixf (region)%acct + dtime
    if (stat_output) then
       statf(region)%acct = statf(region)%acct + dtime
    endif

    ! ----------------------
    ! 3D meteo fields
    ! ----------------------

    ! temperature
    mixf(region)%f3d(temp)%field = mixf(region)%f3d(temp)%field + dtime * temper_dat(region)%data(i1:i2,j1:j2,:) 

    ! specific humidity
    mixf(region)%f3d(hus)%field = mixf(region)%f3d(hus)%field + dtime * humid_dat(region)%data(i1:i2,j1:j2,:) 

    ! air mass (kg -> kg/m2)
    do j = j1, j2
       mixf(region)%f3d(airmass)%field(:,j,:) = mixf(region)%f3d(airmass)%field(:,j,:) + &
            dtime * m_dat(region)%data(i1:i2,j,:) / dxyp(j)
    end do

    ! pressure (already filled above)
    mixf(region)%f3d(pressure)%field = mixf(region)%f3d(pressure)%field + dtime * pres3d

    ! ----------------------
    ! cycle over horizontal domain
    ! ----------------------
    do j = j1, j2
       do i = i1, i2

          ! -----------------------
          ! SURFACE FIELDS
          ! -----------------------
          k = 1

          ! ----------------------
          ! Physical Parameters
          ! ----------------------
          ! surface pressure [Pa]
          mixf(region)%f2d(ps)%field(i,j) = mixf(region)%f2d(ps)%field(i,j) + dtime * sp_dat(region)%data(i,j,k)
          ! precipitation ([m s-1]  -->  [kg m-2 s-1]  with density of water `rol`)
          mixf(region)%f2d(precip)%field(i,j) = mixf(region)%f2d(precip)%field(i,j) + dtime * (cp_dat(region)%data(i,j,k) + lsp_dat(region)%data(i,j,k)) * rol 

          ! density for conversion of aerosol mass densities ( ---> 1/m3 )
          ! (conversion factor 1.E-03 is for  g --> kg)
          dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * Rgas) 

          ! ----------------------
          ! Surface Concentrations in [kg m-3]
          ! ----------------------
          ! POM: (AIS + ACS + COS + AII)
          mixf(region)%f2d(sconcoa)%field(i,j) = mixf(region)%f2d(sconcoa)%field(i,j) + dtime * & 
               dens * (rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii))
          ! BC: (AIS + ACS + COS + AII)
          mixf(region)%f2d(sconcbc)%field(i,j) = mixf(region)%f2d(sconcbc)%field(i,j) + dtime * & 
               dens * (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii))
          ! SO4: (NUS + AIS + ACS + COS)
          mixf(region)%f2d(sconcso4)%field(i,j) = mixf(region)%f2d(sconcso4)%field(i,j) + dtime * & 
               dens * (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos))
          ! Dust:
          mixf(region)%f2d(sconcdust)%field(i,j) = mixf(region)%f2d(sconcdust)%field(i,j) + dtime * & 
               dens * (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi))
          ! Seasalt:
          mixf(region)%f2d(sconcss)%field(i,j) = mixf(region)%f2d(sconcss)%field(i,j) + dtime * & 
               dens * (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos))
          ! NO3:
          mixf(region)%f2d(sconcno3)%field(i,j) = mixf(region)%f2d(sconcno3)%field(i,j) + dtime * &
               dens * rm(i,j,k,iNO3_a)
          ! NH4:
          mixf(region)%f2d(sconcnh4)%field(i,j) = mixf(region)%f2d(sconcnh4)%field(i,j) + dtime * &
               dens * rm(i,j,k,iNH4)
          ! MSA: 
          mixf(region)%f2d(sconcmsa)%field(i,j) = mixf(region)%f2d(sconcmsa)%field(i,j) + dtime * &
               dens * rm(i,j,k,iMSA)

          ! ----------------------
          ! Load in [kg m-2]
          ! ----------------------
          do k = 1, lmr 

             ! POM: (AIS + ACS + COS + AII)
             load_tmp = (rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii))
             mixf(region)%f2d(loadoa)%field(i,j) = mixf(region)%f2d(loadoa)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! BC: (AIS + ACS + COS + AII)
             load_tmp = (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii))
             mixf(region)%f2d(loadbc)%field(i,j) = mixf(region)%f2d(loadbc)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! SO4: (NUS + AIS + ACS + COS)
             load_tmp = (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos))
             mixf(region)%f2d(loadso4)%field(i,j) = mixf(region)%f2d(loadso4)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! Dust: (ACS + COS + ACI + COI)
             load_tmp = (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi))
             mixf(region)%f2d(loaddust)%field(i,j) = mixf(region)%f2d(loaddust)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! Seasalt: (ACS + COS)
             load_tmp = (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos))
             mixf(region)%f2d(loadss)%field(i,j) = mixf(region)%f2d(loadss)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! NO3:
             load_tmp = rm(i,j,k,iNO3_a)
             mixf(region)%f2d(loadno3)%field(i,j) = mixf(region)%f2d(loadno3)%field(i,j) + load_tmp * dtime / dxyp(j)

          end do ! k


       end do ! i 
    end do ! j

    if (stat_output) then
       ! --------------------------------------
       ! stations
       ! --------------------------------------

       if (stat_halo_glob.eq.1) then
          CALL UPDATE_HALO( dgrid(region), rm, mass_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          CALL UPDATE_HALO( dgrid(region), m_dat(region)%data, m_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          CALL UPDATE_HALO( dgrid(region), gph_dat(region)%data, gph_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          do imode = 1, nsol
             CALL UPDATE_HALO( dgrid(region), h2o_mode(region,imode)%d3, h2o_mode(region,imode)%halo, status)
             IF_NOTOK_RETURN(status=1) 
          end do
          CALL UPDATE_HALO( dgrid(region), temper_dat(region)%data(:,:,:), temper_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          CALL UPDATE_HALO( dgrid(region), humid_dat(region)%data(:,:,:), humid_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          ! sp_dat has been updated earlier
          CALL UPDATE_HALO( dgrid(region), cp_dat(region)%data(:,:,1), cp_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          CALL UPDATE_HALO( dgrid(region), lsp_dat(region)%data(:,:,1), lsp_dat(region)%halo, status)
          IF_NOTOK_RETURN(status=1)
          CALL UPDATE_HALO_IBAND ( dgrid(region), dxyp(:), region_dat(region)%halo, status)
       endif

       allocate(volume(i1:ie,j1:je,1:lmr))
       do j = j1, je 
          volume(i1:ie,j,1:lmr) = (gph_dat(region)%data(i1:ie,j,2:lmr+1)-gph_dat(region)%data(i1:ie,j,1:lmr)) * dxyp(j)
       end do 

       do istat = 1, nstations
          if (stat_active(istat)) then
      
             ! masses 
             do i = 1, 18
                statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                      aerocom_do_interp_1d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / & 
                      m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat )
                statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                      aerocom_do_interp_0d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / &
                      m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat )
             end do

             ! number densities
             do i = 19, 19+nmod-1
                statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                      aerocom_do_interp_1d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / &
                      volume(i1:ie,j1:je,1:lmr), region, istat )
                statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                      aerocom_do_interp_0d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / &
                      volume(i1:ie,j1:je,1:lmr), region, istat )
             end do

             ! aerosol water of the wet modes
             do i = 19+nmod, 19+nmod+nsol-1 
                statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                      aerocom_do_interp_1d( h2o_mode(region,statf(region)%f1d(i)%mf%itm5)%d3(i1:ie,j1:je,1:lmr) / &
                      m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat )
                statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                      aerocom_do_interp_0d( h2o_mode(region,statf(region)%f1d(i)%mf%itm5)%d3(i1:ie,j1:je,1:lmr) / &
                      m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat )
             end do

             ! nitrate, ammonium and MSA
             do i= 19+nmod+nsol,19+nmod+nsol+3-1
                statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                      aerocom_do_interp_1d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / &
                      m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat )
                statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                      aerocom_do_interp_0d( rm(i1:ie,j1:je,1:lmr,statf(region)%f1d(i)%mf%itm5) / &
                      m_dat(region)%data(i1:ie,j1:je,1:lmr), region, istat )
             end do

             ! temperature
             i=ntracer_1dstat-nextra_1dstat+1
             statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                   aerocom_do_interp_1d( temper_dat(region)%data(i1:ie,j1:je,:), region, istat )
             statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                   aerocom_do_interp_0d( temper_dat(region)%data(i1:ie,j1:je,:), region, istat )

             ! specific humidity
             i=i+1
             statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                   aerocom_do_interp_1d( humid_dat(region)%data(i1:ie,j1:je,:), region, istat ) 
             statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                   aerocom_do_interp_0d( humid_dat(region)%data(i1:ie,j1:je,:), region, istat )

             ! pressure
             i=i+1 
             statf(region)%f1d(i)%field(istat,:) = statf(region)%f1d(i)%field(istat,:) + dtime * &
                   aerocom_do_interp_1d( pres3d(i1:ie,j1:je,:), region, istat )
             statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                   aerocom_do_interp_0d( pres3d(i1:ie,j1:je,:), region, istat )

             ! surface pressure
             i=i+1
             statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                   aerocom_do_interp_surf( sp_dat(region)%data(i1:ie,j1:je,1), region, istat )

             ! total precipitation
             i=i+1
             statf(region)%f0d(i)%field(istat) = statf(region)%f0d(i)%field(istat) + dtime * &
                   aerocom_do_interp_surf( ( cp_dat(region)%data(i1:ie,j1:je,1) + &
                                            lsp_dat(region)%data(i1:ie,j1:je,1) ) *rol, region, istat )

          endif
       end do 

       deallocate(volume)
       
    endif

    ! ---------------------
    ! DO OPTICS IN PARALLEL 
    ! ---------------------
    ! steps needed for that: 
    ! 1) according to the parallelisation there is different container sizes for 
    !    the results of the optics; these have to be allocated correctly 
    !    (aop_output_ext/a/g/add) 
    ! 2) the optics code is called (init/calculate_aop/done), see documentation
    !    in the optics module
    ! 3) the distributed fields (levels/tracers) are reshaped to global fields 
    !    (opt_ext/a/g/add), according to parallelisation 
    ! 4) fields are communicated such that the result contains the right 
    !    total extinctions, albedos, asymmetry factors etc. 
    ! 5) post-calculations (AOD etc.) are done and results are written
    !    to mixf%../statf% containers as well as to the AOD dump files
    ! 6) done... 
    ! ------------ REMARKS
    ! IMPOSSIBLE / TOO EXPENSIVE (from the AEROCOM list of parameters for QUICKLOOK)
    ! - gf @ 90% RH 
    ! - "AOD@550nm SOA", "AOD@550nm BB"

    ! ---------------------------------
    ! fill current M7 state into an appropriate container
    ! and allocate output fields (ext,a,g)

    ! ---------------------------------
    ! ----- A T T E N T I O N ---------
    ! in case for a 'split', we need a field big enough to contain
    ! various contributions 
    ! (to be synchronously changed with optics_aop_calculate!!) 
    ! ncontr is here number of contributors

    ! Total, SO4, BC, OC, SS, DU, NO3, Water, Fine, Fine Dust, Fine SS
    !ncontr = 10
    ncontr = 11

    ! Additional fields for surface dry aerosol 
    ! to be used for validation against in-situ data: 
    ! extinction, absorption, asymmetry factor
    ! fine-mode contribution to extinction, absorption
    nadd   = 5

    lvec = imr*jmr*lmr

    ! allocate output fields for optics_calculate_aop
    Allocate( aop_output_ext(lvec, size(wdep_out), ncontr ) ) ! extinction
    Allocate( aop_output_a  (lvec, size(wdep_out)         ) ) ! single scattering albedo
    Allocate( aop_output_g  (lvec, size(wdep_out)         ) ) ! asymmetry factor
    Allocate( aop_output_add(lvec, size(wdep_out), nadd   ) ) ! extinction/absorption due to dry aerosol at surface

    call optics_aop_get(lvec, region, size(wdep_out), wdep_out, ncontr, aop_output_ext, aop_output_a, aop_output_g, &
                   status, aop_output_add )
    IF_NOTOK_RETURN(status=1)

    ! allocate here result arrays for ext, abs, ssa, asymmetry parameter, additional parameters (PM10)
    ! ncontr is number of contributors for 'split'
    allocate( opt_ext( i1:i2, j1:j2, lmr, size(wdep_out), ncontr ) ) ; opt_ext = 0.0 
    allocate( opt_a  ( i1:i2, j1:j2, lmr, size(wdep_out)         ) ) ; opt_a   = 0.0 
    allocate( opt_g  ( i1:i2, j1:j2, lmr, size(wdep_out)         ) ) ; opt_g   = 0.0 
    if (stat_output) then
       ! halo cells are needed for interpolation of opt_add
       allocate( opt_add( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lmr, size(wdep_out), nadd   ) ) ; opt_add = 0.0 
    else
       allocate( opt_add( i1:i2, j1:j2, lmr, size(wdep_out), nadd   ) ) ; opt_add = 0.0
    endif

    ! ---------------------------------
    ! unpack results from calculate_aop   
    ! ---------------------------------

    do lwl = 1, size(wdep_out) 
       if( wdep_out(lwl)%split ) then 
          ! fill the array for the extinction coefficients of contributors
          do lct = 1, ncontr 
             opt_ext(:,:,:,lwl,lct) = reshape( aop_output_ext(:,lwl,lct), (/imr,jmr,lmr/) )
          end do
       else
          ! fill only total extinction coefficients
          opt_ext(:,:,:,lwl,1) = reshape( aop_output_ext(:,lwl,1), (/imr,jmr,lmr/) )
       endif 
       opt_a(:,:,:,lwl) = reshape( aop_output_a(:,lwl),(/imr,jmr,lmr/) )
       opt_g(:,:,:,lwl) = reshape( aop_output_g(:,lwl),(/imr,jmr,lmr/) )
       if ( wdep_out(lwl)%insitu ) then
          ! additional fields for split
          do lct = 1, nadd
             opt_add(i1:i2,j1:j2,:,lwl,lct) = reshape( aop_output_add(:,lwl,lct),(/imr,jmr,lmr/) )
          end do
       endif
    end do ! lwl

    ! free temporary arrays for results from calculate_aop
    deallocate( aop_output_ext ) 
    deallocate( aop_output_a  )
    deallocate( aop_output_g  ) 
    deallocate( aop_output_add ) 

    ! the following sections assume that for 550nm there is splitted information 
    ! available and that there is *NO* split for other wavelengths! 
    if( count( (wdep_out(:)%wl == 0.550) .and. wdep_out(:)%split ) /= 1 ) then
       write(gol,*) 'user_output_aerocom: wrong setup with splitted AOD information.'; call goErr
       TRACEBACK
       status=1; return
    end if

    ! -------------------------------------
    ! here additional calculations and 
    ! file writing begin
    ! -------------------------------------

    ! cycle over wavelengths
    do lwl = 1, size(wdep_out) 

       if( wdep_out(lwl)%split ) then 
          ! for 550nm: 
          ! 1) AODs (35 - 44)
          !    fill for contributors (total, SO4, BC, POM, SS, DU, NO3, H2O, fine, fine dust, fine SS)
          ! 2) Absorption for 550nm (45)
          !    Extinction is here the sum of scattering and absorption and 
          !    the scattering part is given by the albedo (SSA). Thus the 
          !    remainder is absorption: Abs = Ext * (1 - SSA)
          indoffset = od550aer
          allocate(temp2d(i1:i2,j1:j2))
          do lct = 1, ncontr
             temp2d = 0.0
             do j = j1, j2
                do i = i1, i2
                   ! multiply with height elements and sum up
                   tabs = 0.0
                   do l = 1, lmr 
                      temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,lct) * & 
                           (gph_dat(region)%data(i,j,l+1) - & 
                            gph_dat(region)%data(i,j,l  ))
                      if( lct == 1 ) then 
                         ! Absorption: do this only once for the totals
                         tabs = tabs + opt_ext(i,j,l,lwl,lct) * (1.0 - opt_a(i,j,l,lwl)) * &  
                              (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l))
                      end if
                   end do
                   ! Absorption: do this only once for the totals
                   if( lct == 1 ) then 
                      mixf(region)%f2d(abs550aer)%field(i,j) = &
                      mixf(region)%f2d(abs550aer)%field(i,j) + tabs * dtime 
                   end if
                end do
             end do
             ! this here is AODs for different contributors, splitted by volumes
             mixf(region)%f2d(indoffset+lct-1)%field = &
             mixf(region)%f2d(indoffset+lct-1)%field + temp2d * dtime 

          end do
          deallocate( temp2d )

          ! Asymmetry Parameter, weighted by AOD and single scattering albedo at each layer
          allocate(temp2d(i1:i2,j1:j2)) ; temp2d = 0.0
          do j = j1, j2
             do i = i1, i2
                wght = 0.0
                do l = 1, lmr
                   dwght = opt_ext(i,j,l,lwl,1) * (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l  )) * opt_a(i,j,l,lwl)
                   temp2d(i,j) = temp2d(i,j) + opt_g(i,j,l,lwl) * dwght
                   wght = wght + dwght
                end do
                temp2d(i,j) = temp2d(i,j) / wght
             end do
          end do
          mixf(region)%f2d(asyaer)%field = mixf(region)%f2d(asyaer)%field + temp2d * dtime
          deallocate(temp2d)
          
          ! Surface Ambient Aerosol Extinction
          mixf(region)%f2d(ec550aer)%field = &
          mixf(region)%f2d(ec550aer)%field + opt_ext(:,:,1,lwl,1) * dtime

       else
          ! for 440/870/350 nm: 
          ! 1) fill total AOD and AAOD and write to containers
          ! 2) dump AOD fields

          if ( wdep_out(lwl)%wl == 0.440 ) then
             indoffset = od440aer
             indabs = abs440aer
          elseif ( wdep_out(lwl)%wl == 0.550 ) then
             indoffset = od550aer
             indabs = abs550aer
          elseif ( wdep_out(lwl)%wl == 0.870 ) then
             indoffset = od870aer
             indabs = abs870aer
          elseif ( wdep_out(lwl)%wl == 0.350 ) then
             indoffset = od350aer
             indabs = abs350aer
          else
             write(gol,*) 'user_output_aerocom: wrong wavelength given for aerocom diagnostics' ; call goErr
             TRACEBACK
             status=1; return
          end if

          ! fill AODs
          allocate(temp2d(i1:i2,j1:j2))
          temp2d = 0.0
          do j = j1, j2
             do i = i1, i2
                ! multiply with height elements and sum up
                tabs = 0.0
                do l = 1, lmr
                   temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,1) * &
                                 (gph_dat(region)%data(i,j,l+1) - &
                                  gph_dat(region)%data(i,j,l  ))
                   tabs = tabs + opt_ext(i,j,l,lwl,1) * (1.0 - opt_a(i,j,l,lwl)) * &
                              (gph_dat(region)%data(i,j,l+1) - gph_dat(region)%data(i,j,l))
                end do
                mixf(region)%f2d(indabs)%field(i,j) = &
                      mixf(region)%f2d(indabs)%field(i,j) + tabs * dtime 
             end do
          end do

          ! write to container 
          mixf(region)%f2d(indoffset)%field = &
          mixf(region)%f2d(indoffset)%field + temp2d * dtime

          deallocate(temp2d)

       endif

       ! 3-D extinction and absorption (m-1)
       if ( aai_output ) then 
          if ( wdep_out(lwl)%wl == 0.550 .or. wdep_out(lwl)%wl == 0.350 ) then
            if ( wdep_out(lwl)%wl == 0.550 ) then
               indoffset = ec5503Daer
            elseif ( wdep_out(lwl)%wl == 0.350 ) then
               indoffset = ec3503Daer
            endif
            mixf(region)%f3d(indoffset)%field = &
            mixf(region)%f3d(indoffset)%field + opt_ext(:,:,:,lwl,1) * dtime
            mixf(region)%f3d(indoffset+1)%field = &
            mixf(region)%f3d(indoffset+1)%field + opt_ext(:,:,:,lwl,1) * (1.0 - opt_a(:,:,:,lwl)) * dtime
          endif
       endif

       if ( wdep_out(lwl)%insitu ) then

          if ( wdep_out(lwl)%wl == 0.440 ) then
             indoffset = ec440dryaer
             indoffset_stat = ec440dryaer_stat
             nadd_max=2
          elseif ( wdep_out(lwl)%wl == 0.550 ) then
             indoffset = ec550dryaer
             indoffset_stat = ec550dryaer_stat
             nadd_max=5
          elseif ( wdep_out(lwl)%wl == 0.870 ) then
             indoffset = ec870dryaer
             indoffset_stat = ec870dryaer_stat
             nadd_max=2
          else
             write(gol,*) 'user_output_aerocom: wrong wavelength given for aerocom diagnostics' ; call goErr
             TRACEBACK
             status=1; return
          end if 

          ! Surface dry aerosol extinction & absorption
          do lct = 1, nadd_max
            mixf(region)%f2d(indoffset+lct-1)%field = mixf(region)%f2d(indoffset+lct-1)%field + &
                                                      opt_add(i1:i2,j1:j2,1,lwl,lct) * dtime
          enddo

          if (stat_output) then
             ! --------------------------------------
             ! stations
             ! --------------------------------------
             if (stat_halo_glob.eq.1) then
                do lct = 1, nadd_max             
                   CALL UPDATE_HALO( dgrid(region), opt_add(:,:,:,lwl,lct), l_halo, status)
                   IF_NOTOK_RETURN(status=1) 
                enddo
             endif 

             do istat = 1, nstations
                if (stat_active(istat)) then
                   do lct = 1, nadd_max
                     statf(region)%f1d(indoffset_stat+lct-1)%field(istat,:) = &
                     statf(region)%f1d(indoffset_stat+lct-1)%field(istat,:) + &
                           aerocom_do_interp_1d( opt_add(i1:ie,j1:je,:,lwl,lct), region, istat ) * dtime
                     statf(region)%f0d(indoffset_stat+lct-1)%field(istat) = &
                     statf(region)%f0d(indoffset_stat+lct-1)%field(istat) + &
                           aerocom_do_interp_0d( opt_add(i1:ie,j1:je,:,lwl,lct), region, istat ) * dtime
                   enddo
                endif
             end do
          endif
 
       endif 

    end do

    ! done
    deallocate( opt_ext, opt_a, opt_g, opt_add )
    nullify(rm)
    nullify(dxyp)
    deallocate( pres3d )

    call goLabel() ; status = 0

  end subroutine output_aerocom_step
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    COLLECT_BUDGETS
  !
  ! !DESCRIPTION: This routine does book keeping of the budget values in 
  !               order to extract daily contributions to 
  !               emissions / dry deposition / wet deposition. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine collect_budgets(region, status)
    !
    ! !USES:
    !
#ifndef without_chemistry
    use ebischeme,      only : buddry_dat => buddep_dat
#endif
    use wet_deposition, only : buddep_dat
    use emission_data,  only : budemi_dat
    use budget_global,  only : nbud_vg
    use global_data, only : region_dat
    use chem_param,  only : iducoi, iduacs, iducos, iduaci, isscos, issacs
    use chem_param,  only : ibccos, ibcaii, ibcacs, ibcais, ino3_a, xmair
    use chem_param,  only : iso4nus, iso4acs, iso4cos, iso4ais
    use chem_param,  only : ipomcos, ipomaii, ipomacs, ipomais
    use chem_param,  only : idms, iso2, ntracet, ntrace
    use chem_param,  only : xmso2, xmso4, xmdms, xmpom, xmbc, xmdust, xmnacl
#ifndef without_sedimentation
    use sedimentation, only : buddep_m7_dat
    use sedimentation, only : nsed, ised 
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   29 Nov 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC


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

    real, dimension(:,:,:,:), allocatable :: collect4d
    real, dimension(:,:,:),   allocatable :: collect3d
    real, dimension(:,:),     allocatable :: tmpbud2d
    real, dimension(:),       pointer     :: dxyp
    integer                               :: imr, jmr, lmr, j
    integer                               :: i1, i2, j1, j2
    !>>> TvN
    integer                               :: l, n
    !<<< TvN

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

    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

    ! this is area element 
    dxyp => region_dat(region)%dxyp

    ! --------------
    ! EMISSIONS
    ! --------------

    ! collect emission budget
    allocate( collect4d(i1:i2,j1:j2,nbud_vg,ntracet) ); collect4d = 0.0 

    ! emissions are in [mole]
    ! convert first to kilomole per area [k-mole/m2]
#ifndef without_emission
    do j = j1, j2
       collect4d(:,j,:,:) = budemi_dat(region)%emi(i1:i2,j,:,:) / dxyp(j) * 1.E-03
    end do
#endif
    allocate( tmpbud2d(i1:i2,j1:j2) )

    ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole]

    ! POM (AIS + ACS + COS + AII)
    tmpbud2d = ( sum(collect4d(:,:,:,iPOMais),3) + sum(collect4d(:,:,:,iPOMacs),3) +   &
                 sum(collect4d(:,:,:,iPOMcos),3) + sum(collect4d(:,:,:,iPOMaii),3) ) * xmpom
    mixf(region)%f2d(emioa)%field = tmpbud2d - emission(region)%f2dslast(:,:,1)
    emission(region)%f2dslast(:,:,1) = tmpbud2d

    ! BC (AIS + ACS + COS + AII)
    tmpbud2d = ( sum(collect4d(:,:,:,iBCais ),3) + sum(collect4d(:,:,:,iBCacs ),3) + & 
                 sum(collect4d(:,:,:,iBCcos ),3) + sum(collect4d(:,:,:,iBCaii ),3) ) * xmbc  
    mixf(region)%f2d(emibc)%field = tmpbud2d - emission(region)%f2dslast(:,:,2)
    emission(region)%f2dslast(:,:,2) = tmpbud2d

    ! SO2
    tmpbud2d = sum(collect4d(:,:,:,iSO2   ),3) * xmso2
    mixf(region)%f2d(emiso2)%field = tmpbud2d - emission(region)%f2dslast(:,:,3)
    emission(region)%f2dslast(:,:,3) = tmpbud2d

    ! SO4 (NUS + AIS + ACS + COS)
    tmpbud2d = ( sum(collect4d(:,:,:,iSO4nus),3) + sum(collect4d(:,:,:,iSO4ais),3) + &  
                 sum(collect4d(:,:,:,iSO4acs),3) + sum(collect4d(:,:,:,iSO4cos),3) ) * xmso4  
    mixf(region)%f2d(emiso4)%field = tmpbud2d - emission(region)%f2dslast(:,:,4)
    emission(region)%f2dslast(:,:,4) = tmpbud2d

    ! Dust (ACS + COS + ACI + COI)
    tmpbud2d = ( sum(collect4d(:,:,:,iDUacs ),3) + sum(collect4d(:,:,:,iDUcos ),3) + &  
                 sum(collect4d(:,:,:,iDUaci ),3) + sum(collect4d(:,:,:,iDUcoi ),3) ) * xmdust       
    mixf(region)%f2d(emidust)%field = tmpbud2d - emission(region)%f2dslast(:,:,5)
    emission(region)%f2dslast(:,:,5) = tmpbud2d

    ! DMS
    tmpbud2d = sum(collect4d(:,:,:,iDMS   ),3) * xmdms
    mixf(region)%f2d(emidms)%field = tmpbud2d - emission(region)%f2dslast(:,:,6)
    emission(region)%f2dslast(:,:,6) = tmpbud2d

    ! Seasalt: (ACS + COS)
    tmpbud2d = ( sum(collect4d(:,:,:,iSSacs ),3) + sum(collect4d(:,:,:,iSScos ),3) ) * xmnacl
    mixf(region)%f2d(emiss)%field = tmpbud2d - emission(region)%f2dslast(:,:,7)
    emission(region)%f2dslast(:,:,7) = tmpbud2d

    deallocate( tmpbud2d )
    deallocate( collect4d )

    ! --------------
    ! DRY DEPOSITION
    ! --------------

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

    ! deposition is in [mole]
    ! convert first to kilomole per area [k-mole/m2]
    do j = j1, j2
#ifndef without_chemistry       
       collect3d(:,j,:) = buddry_dat(region)%dry(i1:i2,j,:) / dxyp(j) * 1.E-3
#endif
       ! Add sedimentation at the surface to dry depostion
       ! Sedimentation budgets have to be summed in the vertical
       ! to get the surface contribution.
#ifndef without_sedimentation
       do l = 1, nbud_vg
        do n = 1, nsed
          collect3d(:,j,ised(n)) = collect3d(:,j,ised(n)) + &
            buddep_m7_dat(region)%sed(i1:i2,j,l,n) / dxyp(j) * 1.E-3
        end do
       end do
#endif
    end do


    allocate( tmpbud2d(i1:i2,j1:j2) )

    ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole]

    ! SO2
    tmpbud2d = collect3d(:,:,iSO2) * xmso2
    mixf(region)%f2d(dryso2)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,1)
    drydepos(region)%f2dslast(:,:,1) = tmpbud2d

    ! POM (AIS + ACS + COS + AII)
    tmpbud2d = ( collect3d(:,:,iPOMais ) + collect3d(:,:,iPOMacs ) + & 
                 collect3d(:,:,iPOMcos ) + collect3d(:,:,iPOMaii ) ) * xmpom
    mixf(region)%f2d(dryoa)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,2)
    drydepos(region)%f2dslast(:,:,2) = tmpbud2d

    ! BC (AIS + ACS + COS + AII)
    tmpbud2d = ( collect3d(:,:,iBCais ) + collect3d(:,:,iBCacs ) + & 
                 collect3d(:,:,iBCcos ) + collect3d(:,:,iBCaii ) ) * xmbc
    mixf(region)%f2d(drybc)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,3)
    drydepos(region)%f2dslast(:,:,3) = tmpbud2d

    ! SO4 (NUS + AIS + ACS + COS)
    tmpbud2d = ( collect3d(:,:,iSO4nus) + collect3d(:,:,iSO4ais) + &
                 collect3d(:,:,iSO4acs) + collect3d(:,:,iSO4cos) ) * xmso4
    mixf(region)%f2d(dryso4)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,4)
    drydepos(region)%f2dslast(:,:,4) = tmpbud2d

    ! Dust (ACS + COS + ACI + COI)
    tmpbud2d = ( collect3d(:,:,iDUacs) + collect3d(:,:,iDUcos) + &
                 collect3d(:,:,iDUaci) + collect3d(:,:,iDUcoi) ) * xmdust
    mixf(region)%f2d(drydust)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,5)
    drydepos(region)%f2dslast(:,:,5) = tmpbud2d

    ! DMS
    ! DMS is NOT deposited in TM5 --> the fields will contain zeros 
    tmpbud2d = collect3d(:,:,iDMS) * xmdms
    mixf(region)%f2d(drydms)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,6)
    drydepos(region)%f2dslast(:,:,6) = tmpbud2d

    ! Seasalt: (ACS + COS)
    tmpbud2d = ( collect3d(:,:,iSSacs) + collect3d(:,:,iSScos) ) * xmnacl
    mixf(region)%f2d(dryss)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,7)
    drydepos(region)%f2dslast(:,:,7) = tmpbud2d

    deallocate( tmpbud2d )
    deallocate( collect3d )

    ! --------------
    ! WET DEPOSITION
    ! --------------

    allocate( collect4d (i1:i2,j1:j2,nbud_vg,ntracet) ); collect4d  = 0.0 

    ! deposition is in [mole]
    ! convert first to kilomole per area [k-mole/m2]
    do j = j1, j2
       collect4d(:,j,:,:) = ( buddep_dat(region)%lsp(i1:i2,j,:,:) + buddep_dat(region)%cp(i1:i2,j,:,:) ) / dxyp(j) * 1.E-3
    end do

    allocate( tmpbud2d(i1:i2,j1:j2) )

    ! on the way: convert from kilomole/m2 to kg/m2 via molar mass [g/mole]

    ! POM (AIS + ACS + COS + AII)
    tmpbud2d = ( sum(collect4d(:,:,:,iPOMais),3) + sum(collect4d(:,:,:,iPOMacs),3) + &
                 sum(collect4d(:,:,:,iPOMcos),3) + sum(collect4d(:,:,:,iPOMaii),3) ) * xmpom
    mixf(region)%f2d(wetoa)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,1)
    wetdepos(region)%f2dslast(:,:,1) = tmpbud2d

    ! BC (AIS + ACS + COS + AII)
    tmpbud2d = ( sum(collect4d(:,:,:,iBCais ),3) + sum(collect4d(:,:,:,iBCacs ),3) + & 
                 sum(collect4d(:,:,:,iBCcos ),3) + sum(collect4d(:,:,:,iBCaii ),3) ) * xmbc
    mixf(region)%f2d(wetbc)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,2)
    wetdepos(region)%f2dslast(:,:,2) = tmpbud2d

    ! SO2
    tmpbud2d = sum(collect4d(:,:,:,iSO2   ),3) * xmso2
    mixf(region)%f2d(wetso2)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,3)
    wetdepos(region)%f2dslast(:,:,3) = tmpbud2d

    ! SO4 (NUS + AIS + ACS + COS)
    tmpbud2d = ( sum(collect4d(:,:,:,iSO4nus),3) + sum(collect4d(:,:,:,iSO4ais),3) + &
                 sum(collect4d(:,:,:,iSO4acs),3) + sum(collect4d(:,:,:,iSO4cos),3) ) * xmso4
    mixf(region)%f2d(wetso4)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,4)
    wetdepos(region)%f2dslast(:,:,4) = tmpbud2d

    ! Dust (ACS + COS + ACI + COI)
    tmpbud2d = ( sum(collect4d(:,:,:,iDUacs ),3) + sum(collect4d(:,:,:,iDUcos ),3) + &
                 sum(collect4d(:,:,:,iDUaci ),3) + sum(collect4d(:,:,:,iDUcoi ),3) ) * xmdust 
    mixf(region)%f2d(wetdust)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,5)
    wetdepos(region)%f2dslast(:,:,5) = tmpbud2d

    ! DMS
    ! DMS is NOT deposited in TM5 --> the fields will contain zeros
    tmpbud2d = sum(collect4d(:,:,:,iDMS   ),3) * xmdms
    mixf(region)%f2d(wetdms)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,6)
    wetdepos(region)%f2dslast(:,:,6) = tmpbud2d

    ! Seasalt: (ACS + COS)
    tmpbud2d = ( sum(collect4d(:,:,:,iSSacs ),3) + sum(collect4d(:,:,:,iSScos ),3) ) * xmnacl
    mixf(region)%f2d(wetss)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,7)
    wetdepos(region)%f2dslast(:,:,7) = tmpbud2d

    deallocate( tmpbud2d )
    deallocate( collect4d )

    nullify(dxyp)

    call goLabel() ; status = 0

  end subroutine collect_budgets
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    AEROCOM_GET_STATINTP
  !
  ! !DESCRIPTION: This routine provides interpolation coefficients to 
  !               (yet another) interpolation routine for the stations. 
  !               The interpolation routine is below ( aerocom_do_interp ). 
  !\\
  ! !INTERFACE:
  !
  subroutine aerocom_get_statintp( region, status )
    !
    ! !USES:
    !
    use dims,      only : dy, dx, yref, xref, xbeg, xend, ybeg, yend, xcyc
    use dims,      only : im, jm, lm 
    use Meteodata, only : gph_dat
    use partools,  only : par_reduce
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region 
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status 
    !
    ! !REVISION HISTORY: 
    !    29 Nov 2010 - Achim Strunk - Creation
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real                          :: rt, sheight, dxr, dyr
    real, dimension(0:lm(region)) :: gheight
    integer                       :: i, l 
    integer                       :: i1, i2, j1, j2
 
    ! --- begin -------------------------------

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

    dyr = dy/yref(region)
    dxr = dx/xref(region)

    do i = 1, nstations 
       ! i-coordinates and weights
       rt = (statlist(i)%statlon - float(xbeg(region)))/dxr + 0.5
       sintp(region)%ixl(i)  = int(rt)
       sintp(region)%wixr(i) = rt - sintp(region)%ixl(i)
       sintp(region)%wixl(i) = 1.0 - sintp(region)%wixr(i)
       if ( sintp(region)%ixl(i) == 0 ) then
            sintp(region)%ixl(i) = im(region)
       endif
       sintp(region)%ixr(i)  = sintp(region)%ixl(i) + 1

       ! j-coordinates and weights
       rt = (statlist(i)%statlat - float(ybeg(region)))/dyr + 0.5
       sintp(region)%jyl(i)  = int(rt)
       sintp(region)%jyr(i)  = sintp(region)%jyl(i) + 1
       sintp(region)%wjyr(i) = rt  - sintp(region)%jyl(i)
       sintp(region)%wjyl(i) = 1.0 - sintp(region)%wjyr(i)
       if ( sintp(region)%jyl(i) == 0            ) sintp(region)%jyl(i) = 1
       if ( sintp(region)%jyr(i) == jm(region)+1 ) sintp(region)%jyr(i) = jm(region)

       ! check if left/lower cell of interpolation calculation is in the processor domain
       ! ensures that each station is handled by a single processor 
       if ( (sintp(region)%ixl(i).ge.i1).and.(sintp(region)%ixl(i).le.i2).and. &
            (sintp(region)%jyl(i).ge.j1).and.(sintp(region)%jyl(i).le.j2) ) then

          stat_active(i)=.true.

          if ( (sintp(region)%ixr(i).eq.i2+1).or.(sintp(region)%jyr(i).eq.j2+1) ) then
             stat_halo_loc=1
          endif

       endif
    enddo

    call PAR_REDUCE( stat_halo_loc, 'max', stat_halo_glob, status, all=.true.)
    if (stat_halo_glob.eq.1) then
       CALL UPDATE_HALO( dgrid(region), gph_dat(region)%data, gph_dat(region)%halo, status)
       IF_NOTOK_RETURN(status=1)
    endif

    do i = 1, nstations
       if ( stat_active(i) ) then
          ! --------
          ! interpolate gph to station, use half level (center of box) & surface
          ! gph is from 1 -> lm(region)+1
          ! --------
          gheight(0) = sintp(region)%wixl(i) * sintp(region)%wjyl(i) * gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyl(i),1) + & 
               sintp(region)%wixr(i) * sintp(region)%wjyl(i) * gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyl(i),1) + & 
               sintp(region)%wixr(i) * sintp(region)%wjyr(i) * gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyr(i),1) + & 
               sintp(region)%wixl(i) * sintp(region)%wjyr(i) * gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyr(i),1)
          do l = 1, lm(region)
             gheight(l) = sintp(region)%wixl(i) * sintp(region)%wjyl(i) * sum(gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyl(i),l:l+1))/2. + & 
                  sintp(region)%wixr(i) * sintp(region)%wjyl(i) * sum(gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyl(i),l:l+1))/2. + & 
                  sintp(region)%wixr(i) * sintp(region)%wjyr(i) * sum(gph_dat(region)%data(sintp(region)%ixr(i),sintp(region)%jyr(i),l:l+1))/2. + & 
                  sintp(region)%wixl(i) * sintp(region)%wjyr(i) * sum(gph_dat(region)%data(sintp(region)%ixl(i),sintp(region)%jyr(i),l:l+1))/2.
          end do

          sintp(region)%tght(i) = gheight(0)

          ! --------
          ! determine layer of station
          ! --------
          ! copy temporary the height of station
          sheight = statlist(i)%statalt

          ! do not allow sampleheight below model surface
          ! -> force station height to model surface
          sheight = max(sheight, gheight(0))

          do l = 1, lm(region)
             if( gheight(l) > sheight ) exit
          end do

          ! success??
          if( l > lm(region) ) then 
             ! not successful, station higher than model (strange!!!) 
             write(gol,'("Error in retrieving vertical layer of station")'); call goErr 
             TRACEBACK
             status=1
             return
          endif

          ! restrict lower index to surface
          sintp(region)%lzl(i)  = max(l-1,1)
          sintp(region)%lzr(i)  = l
          sintp(region)%wlzl(i) = (gheight(l) - sheight) / (gheight(l) - gheight(l-1))
          sintp(region)%wlzr(i) = 1.0 - sintp(region)%wlzl(i)

       endif

    end do

    call goLabel() ; status = 0

  end subroutine aerocom_get_statintp
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    AEROCOM_DO_INTERP_1D
  !
  ! !DESCRIPTION: Interpolation of a given 3d field, using module wide
  !               interpolation coeffients and indices for stations (sintp). 
  !\\
  !\\
  ! !INTERFACE:
  !
  function aerocom_do_interp_1d( f3d, region, istat ) 
    !
    ! !INPUT PARAMETERS:
    !
    real, dimension(:,:,:), intent(in)     :: f3d
    integer,                intent(in)     :: region
    integer,                intent(in)     :: istat
    !
    ! !RETURN VALUE:
    !
    real,dimension (levi%nlev) :: aerocom_do_interp_1d
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                     :: i 
    integer                     :: i1, i2, j1, j2
    integer                     :: il,ir,jl,jr

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

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

    i = istat

    ! shift index in the input array 
    il=sintp(region)%ixl(i) - i1 + 1
    ir=sintp(region)%ixr(i) - i1 + 1
    jl=sintp(region)%jyl(i) - j1 + 1
    jr=sintp(region)%jyr(i) - j1 + 1

       ! ordering 
       ! i = llllrrrr, j = llrrllrr
       aerocom_do_interp_1d(:) = &
            sintp(region)%wixl(i) * sintp(region)%wjyl(i) * f3d(il,jl,:) + & 
            sintp(region)%wixl(i) * sintp(region)%wjyr(i) * f3d(il,jr,:) + & 
            sintp(region)%wixr(i) * sintp(region)%wjyl(i) * f3d(ir,jl,:) + & 
            sintp(region)%wixr(i) * sintp(region)%wjyr(i) * f3d(ir,jr,:)

    return

  end function aerocom_do_interp_1d
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    AEROCOM_DO_INTERP_0D
  !
  ! !DESCRIPTION: Interpolation of a given 3d field, using module wide
  !               interpolation coeffients and indices for stations (sintp). 
  !\\
  !\\
  ! !INTERFACE:
  !
  function aerocom_do_interp_0d( f3d, region, istat )
    !
    ! !INPUT PARAMETERS:
    !
    real, dimension(:,:,:), intent(in)     :: f3d
    integer,                intent(in)     :: region
    integer,                intent(in)     :: istat
    !
    ! !RETURN VALUE:
    !
    real :: aerocom_do_interp_0d
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                     :: i
    integer                     :: i1, i2, j1, j2
    integer                     :: il,ir,jl,jr

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

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

    i = istat

    ! shift index in the input array 
    il=sintp(region)%ixl(i) - i1 + 1
    ir=sintp(region)%ixr(i) - i1 + 1
    jl=sintp(region)%jyl(i) - j1 + 1
    jr=sintp(region)%jyr(i) - j1 + 1

       ! ordering 
       ! i = llllrrrr, j = llrrllrr, l = lrlrlrlr
       aerocom_do_interp_0d = &
            sintp(region)%wixl(i) * sintp(region)%wjyl(i) * sintp(region)%wlzl(i) * f3d(il,jl,sintp(region)%lzl(i)) + &
            sintp(region)%wixl(i) * sintp(region)%wjyl(i) * sintp(region)%wlzr(i) * f3d(il,jl,sintp(region)%lzr(i)) + &
            sintp(region)%wixl(i) * sintp(region)%wjyr(i) * sintp(region)%wlzl(i) * f3d(il,jr,sintp(region)%lzl(i)) + &
            sintp(region)%wixl(i) * sintp(region)%wjyr(i) * sintp(region)%wlzr(i) * f3d(il,jr,sintp(region)%lzr(i)) + &
            sintp(region)%wixr(i) * sintp(region)%wjyl(i) * sintp(region)%wlzl(i) * f3d(ir,jl,sintp(region)%lzl(i)) + &
            sintp(region)%wixr(i) * sintp(region)%wjyl(i) * sintp(region)%wlzr(i) * f3d(ir,jl,sintp(region)%lzr(i)) + &
            sintp(region)%wixr(i) * sintp(region)%wjyr(i) * sintp(region)%wlzl(i) * f3d(ir,jr,sintp(region)%lzl(i)) + &
            sintp(region)%wixr(i) * sintp(region)%wjyr(i) * sintp(region)%wlzr(i) * f3d(ir,jr,sintp(region)%lzr(i))

    return

  end function aerocom_do_interp_0d
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:    AEROCOM_DO_INTERP_SURF
  !
  ! !DESCRIPTION: Interpolation of a given surface field, using module wide
  !               interpolation coeffients and indices for stations (sintp). 
  !\\
  !\\
  ! !INTERFACE:
  !
    function aerocom_do_interp_surf( f2d, region, istat )
    !
    ! !INPUT PARAMETERS:
    !
    real, dimension(:,:),   intent(in)     :: f2d
    integer,                intent(in)     :: region
    integer,                intent(in)     :: istat
    !
    ! !RETURN VALUE:
    !
    real :: aerocom_do_interp_surf
    !
    ! !REVISION HISTORY: 
    !    6 Feb 2011 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                     :: i
    integer                     :: i1, i2, j1, j2
    integer                     :: il,ir,jl,jr

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

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

    i = istat

    ! shift index in the input array 
    il=sintp(region)%ixl(i) - i1 + 1
    ir=sintp(region)%ixr(i) - i1 + 1
    jl=sintp(region)%jyl(i) - j1 + 1
    jr=sintp(region)%jyr(i) - j1 + 1

       ! ordering 
       ! i = llllrrrr, j = llrrllrr
       aerocom_do_interp_surf = &
            sintp(region)%wixl(i) * sintp(region)%wjyl(i) * f2d(il,jl) + &
            sintp(region)%wixl(i) * sintp(region)%wjyr(i) * f2d(il,jr) + &
            sintp(region)%wixr(i) * sintp(region)%wjyl(i) * f2d(ir,jl) + &
            sintp(region)%wixr(i) * sintp(region)%wjyr(i) * f2d(ir,jr)

    return

  end function aerocom_do_interp_surf
  !EOC

end module user_output_aerocom