ecearth3/sources/tm5mp/proj/bacchus/user_output_general.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_GENERAL
!
! !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_GENERAL
  !
  ! !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  
  use mo_aero_m7 , only: ncomp,nmod,naermod
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: output_general_init
  public :: output_general_step
  public :: output_general_write
  public :: output_general_done
  public :: wdep_out

  ! Control variables
  character(len=20), public :: gen_exper ! general experiment name
  character(len=20), public :: gen_freq  ! general output frequency
  logical,public            :: all_chemistry=.false.
  logical, parameter :: aai_output=.false.
  !
  ! diagnostic CCN
  logical, public ::  nCCNdiag   = .false.     ! diagnostic CCN under given supersaturations for comparison to CCNC measurements
                                       !   nCCNdiag = .false. no diagnostic CCN
                                       !   nCCNdiag = .true.  with diagnostic CCN
  integer,public :: nSat = 0                  ! Number of supersaturations for which number of CCN will be calculated
  real, public,  allocatable :: SuperSat(:)

  integer,public :: n_rad = 3                  ! bacchus output at n diam
  integer, public :: HG_rad(3)=(/50,80, 120/)
  integer, public :: N_diam(3)=(/50,80, 120/)
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*),  parameter ::  mname = 'user_output_general'

  ! Changed naming convections to AeroCom Control 2015 experiment
  character(len=20), parameter :: f_dataid='general', f_model='TM5'
  character(len=20), parameter :: f_dimmix='global',   f_dimstat='stations' 
  !
  ! !PRIVATE TYPES:
  !
  type metafields
     integer                       :: itm5
     character(len=20)             :: 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
     logical                         :: writeout
 end type field3d

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

  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
     real,dimension(:),pointer                :: lon             ! x dimension of requested field
     real,dimension(:),pointer                :: lat             ! y dimension of requested field
     real,dimension(:),pointer               :: lev             ! z dimension of requested field
     integer               :: lonid             ! x dimension id in nc
     integer               :: latid             ! y dimension id in nc
     integer               :: levid             ! z dimension id in nc
     integer               :: timeid             ! time dimension id in nc
     integer               :: time_varid
  end type mixfile
  type(mixfile), dimension(nregions), save  :: mixf  
  

  integer::lon_varid
  integer::lat_varid
  integer::lev_varid

  ! diagnostic CCN
  integer, allocatable :: ind_CCN(:)
  integer, allocatable :: ind_HG(:)
  integer, allocatable :: ind_N(:)

  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                         ::  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          ::  time_reftime6(6) = (/2000,01,01,00,00,00/)
  character(len=*), parameter   ::  time_units       = 'hours since 2000-01-01 00:00' 

  ! sizes of arrays
  integer, save :: ntracer_3d, ntracer_2d
  integer, save :: ntracer_1dstat, ntracer_0dstat, nstations 
  integer, save :: nextra_1dstat
  integer, save :: n_2d_vars=0
  integer, save :: n_3d_vars=0
  ! index pointers for 1d/2d/3d fields in mixf
  integer, save :: temp, hus, airmass, pressure
  integer, save :: p_gas_so4,p_liq_so4,p_elvoc,p_svoc,d_nuc,m_nuc_su,m_nuc_soa,p_svoc2d,p_elvoc2d,p_gas_so42d,p_liq_so42d
  integer, save :: p_elohterp,p_elo3terp,p_elohisop,p_elo3isop,p_svohterp,p_svo3terp,p_svohisop,p_svo3isop

  integer, save :: coag1, gr1_2, co1_soa, co1_su, co2_soa, co3_soa, co4_soa, co5_soa,gph
  integer, save :: so4pm1,bcpm1,orgpm1,soapm1,dupm1,sspm1,nh4pm1,no3pm1
  integer, save :: so4,bc,org,du,ss,nh4,no3,n3,mmf1,mmf2,mmf3,mmf4
  integer, save :: ec5503Daer, abs5503Daer, ec3503Daer, abs3503Daer
  integer, save :: ps          , precip      , sconcoa     , sconcsoa    , sconcbc     , sconcso4    , sconcdust   
  integer, save :: sconcss     , sconcno3    , loadoa      , loadsoa     , 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       , wetsoa      , wetbc       , wetso2      
  integer, save :: emiterp     , emiisop
  integer, save :: wetso4      , wetdust     , wetdms      , wetss       , od550aer    , od550so4, od550soa    
  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       , drysoa      , sconcnh4
  integer, save :: abs440aer   , ec440dryaer , abs440dryaer
  integer, save :: abs870aer   , ec870dryaer , abs870dryaer
  integer, save :: od350aer    , abs350aer,ccn02
  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
  integer, save :: tr2d_1, tr2d_2, tr2d_3, tr2d_4, tr2d_5, tr2d_6, tr2d_7, tr2d_8, tr2d_9, tr2d_10
  integer, save :: tr2d_11, tr2d_12, tr2d_13, tr2d_14, tr2d_15, tr2d_16, tr2d_17, tr2d_18, tr2d_19
  integer, save :: tr2d_20, tr2d_21, tr2d_22, tr2d_23
  integer, save :: cc01, cc02, cc03, cc04, cc05, cc06, cc07
  integer, save :: cc3d01, cc3d02, cc3d03, cc3d04, cc3d05, cc3d06, cc3d07
  integer, save :: rw01, rw02, rw03, rw04, rw05, rw06, rw07
  integer, save :: rd01, rd02, rd03, rd04, rd05, rd06, rd07
  integer, save :: h2o1, h2o2, h2o3, h2o4
  integer, save :: tr3d_1, tr3d_2, tr3d_3, tr3d_4, tr3d_5, tr3d_6, tr3d_7, tr3d_8, tr3d_9, tr3d_10
  integer, save :: tr3d_11, tr3d_12, tr3d_13, tr3d_14, tr3d_15, tr3d_16, tr3d_17, tr3d_18, tr3d_19
  integer, save :: rw3d01, rw3d02, rw3d03, rw3d04, rw3d05, rw3d06, rw3d07
  integer, save :: rd3d01, rd3d02, rd3d03, rd3d04, rd3d05, rd3d06, rd3d07
  integer, save :: h2o3d1, h2o3d2, h2o3d3, h2o3d4
  integer, save ,dimension(2*nmod)   :: radius
  integer, save ,dimension(ncomp)    :: load
  integer, save ,dimension(ncomp)    :: drydep
  integer, save ,dimension(ncomp)    :: emi
  integer, save ,dimension(ncomp)    :: wetdep
  integer, save ,dimension(ncomp)    :: sed
  integer, save ,dimension(nmod)     :: number
  integer, save ,dimension(naermod)  :: masses
  integer, save :: ngas_output=10
  integer, save, dimension(:),allocatable :: gas_output

  integer::itim_opt_out
  !
  ! !REVISION HISTORY: 
  !   16 Nov 2010 - Achim Strunk
  !   8  Dec 2016 - Tommi Bergman  - Modified from AEROCOM-output
  !
  ! !REMARKS:
  ! (1) compiled only if with_m7 is used.
  !
  !EOP
  !------------------------------------------------------------------------

contains


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_GENERAL_INIT
  !
  ! !DESCRIPTION: Initialise various parameters, eg.,
  !               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_general_init(status, iregion)
    !
    ! !USES:
    !
    !use MeteoData,     only : sp_dat, set
    use GO,            only: GO_Timer_Def, GO_Timer_End, GO_Timer_Start
    use chem_param,    only : ntracet,ntrace,names,mode_tracers,mode_start,mode_nm
    use chem_param,    only : inus_n, iso4nus, isoanus,nmod
    use chem_param,    only : iais_n, iso4ais, ibcais, ipomais, isoaais 
    use chem_param,    only : iacs_n, iso4acs, ibcacs, ipomacs, isoaacs, issacs, iduacs
    use chem_param,    only : icos_n, iso4cos, ibccos, ipomcos, isoacos, isscos, iducos
    use chem_param,    only : iaii_n,          ibcaii, ipomaii, isoaaii
    use chem_param,    only : iaci_n, icoi_n,                           iduaci, iducoi
    use chem_param,    only : ino3_a, inh4, imsa
    use chem_param,    only : Kap_su, Kap_pom, Kap_soa, Kap_bc, Kap_du, Kap_ss, Kap_na2so4, Kap_no3, Kap_msa
    !use chem_param,    only : io3, idms, imsa, inh4, iterp, ioh, ino3, ielvoc, isvoc, iisop,iso2,iso4
    use chem_param,    only : io3,  ih2o2,  ieth,  irooh, idms,   imsa,   imgly, inox,      ico,     ipar,    ipan    
    use chem_param,    only : ich2o, iald2,  iisop, iso4,  irn222, ino3_a, ich4,  ich3o2h,   iole,    iorgntr, inh3    
    use chem_param,    only : ic2h6, iethoh, ic3h8, iterp, iacet,  iispd,  ihono, ich3o2no2, ino,     iho2,    ich3o2
    use chem_param,    only : in2o5, ihno4,  ic2o3, iror,  irxpar, ixo2,   ixo2n, inh2,      ih2opart,ic3h7o2, ihypro2
    use chem_param,    only : isvoc, iso2,  inh4,   ipb210,io3s,  ich3oh, ihcooh ,  ioh,   ino2,   ino3, iaco2, inh2o2, ielvoc
    USE mo_aero_m7,   ONLY : nmod, nsol
    use dims,          only : nregions, newsrun, idatee, idatei ,sdate_simulation 
    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, idatei

    use NetCDF, only : NF90_Def_Var,NF90_put_var
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          :: status
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in), optional :: iregion 
    !
    ! !REVISION HISTORY: 
    !      Nov 2010 - Achim Strunk
    !   16 Dec 2016 - Tommi Bergman KNMI 
    !         - Modified from AEROCOM output
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    integer               :: imr, jmr, lmr, access_mode_sta
    integer               :: lat_dimid,lon_dimid,lev_dimid,kappa_dimid
    integer               :: year_varid,month_varid,kso4_varid, kss_varid,kpom_varid,kbc_varid,kdu_varid,kno3_varid,kna2so4_varid,ksoa_varid,kmsa_varid
    integer               :: region, varid
    integer               :: io, istat, i, j, n, sc, iSat
    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, time_shift
    real                  :: dlat, dlon
    integer               :: mlength    
    integer               :: itracer,index_aer,imass,iout
    character(len=14)     :: str, str2
    character(len=28)     :: str_reftime
    ! --- 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,'("General 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 

       call GO_Timer_Def( itim_opt_out, 'optics-output', status )
       IF_NOTOK_RETURN(status=1)

       ! 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 = 160+nSat
       ! list removed, not needed anymore.

       ! 2d-output list needed for the moment
       ! optics needs to be changed first, since it relies on having optical  variables 
       ! using indices 38->
       ntracer_2d   = 122
       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
       od550soa     = 42  ! 2d
       od550ss      = 43  ! 2d
       od550dust    = 44  ! 2d
       od550no3     = 45  ! 2d
       od550aerh2o  = 46  ! 2d
       od550lt1aer  = 47  ! 2d
       od550lt1dust = 48  ! 2d
       od550lt1ss   = 49  ! 2d
       ! --- end split order
       abs550aer    = 50  ! 2d
       asyaer       = 51  ! 2d
       ec550aer     = 52  ! 2d
       ! --- begin in-situ data order
       ec550dryaer  = 53  ! 2d
       abs550dryaer = 54  ! 2d
       asydryaer    = 55  ! 2d
       ec550drylt1aer  = 56  ! 2d
       abs550drylt1aer = 57  ! 2d
       ! --- end in-situ data order
       !
       od440aer     = 58  ! 2d
       abs440aer    = 59  ! 2d
       ec440dryaer  = 60  ! 2d
       abs440dryaer = 61  ! 2d
       !
       od870aer     = 62  ! 2d
       abs870aer    = 63  ! 2d
       ec870dryaer  = 64  ! 2d
       abs870dryaer = 65  ! 2d
       !
       od350aer     = 66  ! 2d
       abs350aer    = 67  ! 2d
       !
       tr2d_1=68
       tr2d_2=69
       tr2d_3=70
       tr2d_4=71
       tr2d_5=72
       tr2d_6=73
       tr2d_7=74
       tr2d_8=75
       tr2d_9=76
       tr2d_10=77
       tr2d_11=78
       tr2d_12=79
       tr2d_13=80
       tr2d_14=81
       tr2d_15=82
       tr2d_16=83
       tr2d_17=84
       tr2d_18=85
       tr2d_19=86
       tr2d_20=87
       tr2d_21=88
       tr2d_22=89
       tr2d_23=90
       cc01=91
       cc02=92
       cc03=93
       cc04=94
       cc05=95
       cc06=96
       cc07=97
       h2o1=98
       h2o2=99
       h2o3=100
       h2o4=101
       rw01=102
       rw02=103
       rw03=104
       rw04=105
       rw05=106
       rw06=107
       rw07=108
       rd01=109
       rd02=110
       rd03=111
       rd04=112
       sconcsoa=113
       loadsoa=114
       drysoa=115
       wetsoa=116
       emiterp=117
       emiisop=118
       p_svoc2d=119
       p_elvoc2d=120
       p_liq_so42d=121
       p_gas_so42d=122
    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_general_init: nbud_vg /= lmr'; call goErr
          write(gol,*)'output_general_init: YOU MUST ADD THE PROJ/BUDGET10 TO SRC CODE !!!!!'; call goErr
          TRACEBACK
          status=1; return
       end if

       !  initialize the output:
       if( newsrun ) then 
          if (all_chemistry) then
             ngas_output=57
             allocate(gas_output(ngas_output))
             gas_output=(/&
                  io3,   ih2o2,  ieth,  irooh, idms,   imsa,   imgly, inox,      ico,     ipar,    ipan,    iso2,  inh4,   ipb210, &
                  ich2o, iald2,  iisop, iso4,  irn222, ino3_a, ich4,  ich3o2h,   iole,    iorgntr, inh3,    io3s,  ich3oh, ihcooh, &
                  ic2h6, iethoh, ic3h8, iterp, iacet,  iispd,  ihono, ich3o2no2, ino,     iho2,    ich3o2,  ioh,   ino2,   ino3,   &
                  in2o5, ihno4,  ic2o3, iror,  irxpar, ixo2,   ixo2n, inh2,      ih2opart,ic3h7o2, ihypro2, iaco2, inh2o2, ielvoc, &
                  isvoc/)
          else
             allocate(gas_output(ngas_output))
             gas_output=(/io3,idms,iso4,iso2,iterp,ioh,ino3,ielvoc,isvoc,iisop /)
          end if
          call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
          !timeidx=1
          ! target structures for output
          allocate(mixf (region)%f3d(ntracer_3d))
          allocate(mixf (region)%f2d(ntracer_2d))
          
          ! allocate data holder for maximum amount of tracers
          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

          ! METEO variables
          call add_Variable(region,-1     , 'temp        ', 'Temperature                        ' , 'K         ',3,temp,status  )
          call add_Variable(region,-1     , 'hus         ', 'Specific Humidity                  ' , '1         ',3,hus,status  )
          call add_Variable(region,-1     , 'airmass     ', 'Air Mass                           ' , 'kg m-2    ',3,airmass,status  )
          call add_Variable(region,-1     ,'pres         ', 'Pressure                           ' , 'Pa        ',3,pressure,status  )
          !call add_Variable(region,-1     ,'prod_gas_so4 ', 'production of so4                  ' , 'kg m-2s-1',3,p_gas_so4,status  )
          !call add_Variable(region,-1     ,'prod_liq_so4 ', 'production of so4                  ' , 'kg m-2s-1',3,p_liq_so4,status  )
          !call add_Variable(region,-1     ,'prod_elvoc   ', 'production of elvoc                ' , 'kg m-2s-1',3,p_elvoc,status  )
          !call add_Variable(region,-1     ,'prod_svoc    ', 'production of svoc                 ' , 'kg m-2s-1',3,p_svoc,status  )
         
          !call add_Variable(region,-1     ,'p_el_ohterp  ', 'production of elvoc, ohterp        ' , 'kg m-2s-1',3,p_elohterp,status  )
          !call add_Variable(region,-1     ,'p_el_o3terp  ', 'production of elvoc, o3terp        ' , 'kg m-2s-1',3,p_elo3terp,status  )
          !call add_Variable(region,-1     ,'p_el_ohisop  ', 'production of elvoc, ohisop        ' , 'kg m-2s-1',3,p_elohisop,status  )
          !call add_Variable(region,-1     ,'p_el_o3isop  ', 'production of elvoc, o3isop        ' , 'kg m-2s-1',3,p_elo3isop,status  )
          !call add_Variable(region,-1     ,'p_sv_ohterp  ', 'production of svoc, ohterp         ' , 'kg m-2s-1',3,p_svohterp,status  )
          !call add_Variable(region,-1     ,'p_sv_o3terp  ', 'production of svoc, o3terp         ' , 'kg m-2s-1',3,p_svo3terp,status  )
          !call add_Variable(region,-1     ,'p_sv_ohisop  ', 'production of svoc, ohisop         ' , 'kg m-2s-1',3,p_svohisop,status  )
          !call add_Variable(region,-1     ,'p_sv_o3isop  ', 'production of svoc, o3isop         ' , 'kg m-2s-1',3,p_svo3isop,status  )
          !call add_Variable(region,-1     ,'d_nuc        ', 'production of particles            ' , ' cm-3s-1',3,d_nuc,status  )
          !call add_Variable(region,-1     ,'m_nuc_su     ', 'nucleation sulfate                 ' , ' molec cm-3s-1',3,m_nuc_su,status  )
          !call add_Variable(region,-1     ,'m_nuc_soa    ', 'nucleation ELVOC                   ' , ' µg m-3s-1',3,m_nuc_soa,status  )
          !call add_Variable(region,-1     ,'gr1_2        ', 'growth of mode 1 into 2            ' , ' cm-3s-1',3,gr1_2,status  )
          !call add_Variable(region,-1     ,'coag1        ', 'coagulation sink of mode 1         ' , ' cm-3s-1',3,coag1,status  )
          !call add_Variable(region,-1     ,'co1_soa      ', 'condensation of ELVOC to mode 1    ' , ' µg m-3s-1',3,co1_soa,status  )
          !call add_Variable(region,-1     ,'co2_soa      ', 'condensation of ELVOC to mode 2    ' , ' µg m-3s-1',3,co2_soa,status  )
          !call add_Variable(region,-1     ,'co3_soa      ', 'condensation of ELVOC to mode 3    ' , ' µg m-3s-1',3,co3_soa,status  )
          !call add_Variable(region,-1     ,'co4_soa      ', 'condensation of ELVOC to mode 4    ' , ' µg m-3s-1',3,co4_soa,status  )
          !call add_Variable(region,-1     ,'co5_soa      ', 'condensation of ELVOC to mode 5    ' , ' µg m-3s-1',3,co5_soa,status  )
          !call add_Variable(region,-1     ,'co1_su       ', 'condensation of sulfate to mode 1  ' , ' molec cm-3s-1',3,co1_su,status  )
          
          call add_Variable(region,-1     ,'mf1       ', 'mass of so4 in pm1  ' , ' ug m-3',3,mmf1,status  )
          call add_Variable(region,-1     ,'mf2       ', 'mass of BC in pm1  ' , ' ug m-3',3,mmf2,status  )
          call add_Variable(region,-1     ,'mf3       ', 'mass of ORG in pm1  ' , ' ug m-3',3,mmf3,status  )
          call add_Variable(region,-1     ,'mf4       ', 'mass of SOA in pm1  ' , ' ug m-3',3,mmf4,status  )

          call add_Variable(region,-1     ,'SO4pm1       ', 'mass of so4 in pm1  ' , ' ug m-3',3,so4pm1,status  )
          call add_Variable(region,-1     ,'BCpm1       ', 'mass of BC in pm1  ' , ' ug m-3',3,bcpm1,status  )
          call add_Variable(region,-1     ,'ORGpm1       ', 'mass of ORG in pm1  ' , ' ug m-3',3,orgpm1,status  )
          call add_Variable(region,-1     ,'SOApm1       ', 'mass of SOA in pm1  ' , ' ug m-3',3,soapm1,status  )
          call add_Variable(region,-1     ,'SSpm1       ', 'mass of SS in pm1  ' , ' ug m-3',3,sspm1,status  )
          call add_Variable(region,-1     ,'DUpm1       ', 'mass of DU in pm1  ' , ' ug m-3',3,dupm1,status  )
          
          call add_Variable(region,-1     ,'NH4pm1       ', 'mass of NH4 in pm1  ' , ' ug m-3',3,nh4pm1,status  )
          
          call add_Variable(region,-1     ,'NO3pm1       ', 'mass of NO3 in pm1  ' , ' ug m-3',3,no3pm1,status  )
          
          call add_Variable(region,-1     ,'SO4       ', 'mass of SO4   ' , ' ug S m-3',3,so4,status  )
          call add_Variable(region,-1     ,'BC       ', 'mass of BC   ' , ' ug C m-3',3,bc,status  )
          call add_Variable(region,-1     ,'ORG       ', 'mass of ORG  ' , ' ug C m-3',3,org,status  )
          call add_Variable(region,-1     ,'SS       ', 'mass of SS    ' , ' ug m-3',3,ss,status  )
          call add_Variable(region,-1     ,'DU       ', 'mass of DU    ' , ' ug m-3',3,du,status  )
          
          call add_Variable(region,-1     ,'NH4       ', 'mass of NH4    ' , ' ug m-3',3,nh4,status  )
          
          call add_Variable(region,-1     ,'NO3       ', 'mass of NO3    ' , ' ug m-3',3,no3,status  )
          
          
          IF (nCCNdiag) THEN
             ALLOCATE(ind_CCN(nSat))
             DO iSat=1,nSat
                !write(1111,*)SuperSat(iSat)*1000.e0 
                if (supersat(isat)<0.0006)then
                   !WRITE(str,'("CCN",I0.3)') SuperSat(iSat)*10000.e0 !CCN005 (0.0005*1e4=5)
                   WRITE(str,'("CCN005")') !SuperSat(iSat)*10000.e0 !CCN005 (0.0005*1e4=5)
                else if (supersat(isat)<0.0008)then
                   WRITE(str,'("CCN0075")') !SuperSat(iSat)*100000.e0!CCN0075 (0.00075*1e5=75
                else if (supersat(isat)<0.0012)then
                   WRITE(str,'("CCN01")') !SuperSat(iSat)*1000.e0 !CCN01 (0.001*1e3=1)
                else if (supersat(isat)<0.0016)then
                   WRITE(str,'("CCN015")') !SuperSat(iSat)*10000.e0!CCN015 0.0015*1e4=15)
                else 
                   !write(1112,*)SuperSat(iSat)*1000.e0 
                   WRITE(str,'("CCN",I0.2)') nint(SuperSat(iSat)*1000.e0)  !CCN0[23456789] & CCN10
                end if
                !write(1112,*) str , int(SuperSat(iSat)*1000.e0) , nint(SuperSat(iSat)*1000.e0), (SuperSat(iSat)*1000.e0)
                WRITE(str2,'("CCN at ",F4.2,"%")') SuperSat(iSat)*100.e0
                call add_Variable(region,-1     ,str, str2, ' cm-3',3,ind_CCN(iSat),status  )
             ENDDO
             allocate(ind_HG(n_rad))
             allocate(ind_N(n_rad))
             
             do isat=1,n_rad
                if (HG_rad(isat)>99) then
                   WRITE(str,'("HG",I3)') HG_rad(isat)
                else
                   WRITE(str,'("HG",I2)') HG_rad(isat)
                end if
                WRITE(str2,'("HG at ",I3,"nm")') HG_rad(isat)
                call add_variable(region,-1     ,str, str2, '1',3,ind_HG(iSat),status )
                
                if (N_diam(isat)>99)then
                   WRITE(str,'("N",I3)') N_diam(isat)
                else
                   WRITE(str,'("N",I2)') N_diam(isat)
                end if
                !write(1000,*)str
                WRITE(str2,'("N at ",I3,"nm")') N_diam(isat)
                call add_variable(region,-1     ,str, str2, ' cm-3',3,ind_N(iSat),status )
             end do
          ENDIF
          call add_variable(region,-1     ,'N3', 'N at 3nm', ' cm-3',3,n3,status )
          call add_Variable(region,-1     ,'gph', 'Geopotential height', 'm',3,gph,status  )
          call add_Variable(region,-1     ,'CCN', 'CCN at 0.2%', '#/cm3',3,ccn02,status  )
          ! mass concentrations of different modes
          imass=0
          do i=1,nmod
             do j=1,mode_nm(i)!ncomp
                !get tracer id from mode_tracers array
                itracer=mode_tracers(j,i)
                ! if there is a tracer, add variable to output array
                !if (itracer>0)then
                imass=imass+1
                call add_Variable(region,itracer     ,'M_'//names(itracer), 'mass concentration of '//names(itracer), 'kg m-3',3,index_aer,status  )
                !end if
                ! save output array id for masses
                masses(imass)=index_aer
             end do
          end do
          call add_Variable(region,imsa     ,'M_MSA', 'mass concentration of MSA (ACS)', 'kg m-3',3,index_aer,status  )
          call add_Variable(region,ino3_a    ,'M_NO3', 'mass concentration of NO3 (ACS)', 'kg m-3',3,index_aer,status  )
          call add_Variable(region,inh4    ,'M_NH4', 'mass concentration of NH4 (ACS)', 'kg m-3',3,index_aer,status  )
          
          ! Mode number concentrations
          do i=1,nmod
             !get tracer id from mode_tracers array
             ! 0-index of mode_tracers has the number for each mode
             itracer=mode_tracers(0,i)
             ! if there is a tracer, add variable to output array
             if (itracer>0)then
                ! 
                call add_Variable(region,itracer     ,'N_'//names(itracer)(1:3), 'number concentration of '//names(itracer), '1 m-3',3,index_aer,status  )
                
             end if
             ! save output array id for number
             number(i)=index_aer
          end do
          
          ! gas-phase output
!!$          do i=1,ngas_output
!!$             itracer=gas_output(i)
!!$             call add_Variable(region,itracer     ,'GAS_'//names(itracer), 'gas phase concentration of '//names(itracer), 'kg m-3',3,index_aer,status  )
!!$          end do
!!$
!!$          ! particle water
!!$          do i=1,nsol
!!$             !write(iout,'(I2.2)') i
!!$             itracer=mode_start(i)
!!$             call add_Variable(region,-1     , 'aerh2o3d_'//names(itracer)(1:3), 'concentration of aerosol water of mode '//names(itracer)(1:3) ,  'kg m-3',3,index_aer,status  )
!!$
!!$          end do

          ! Radius of the particles
          ! wet radius for all
!!$          do i=1,nmod
!!$             itracer=mode_start(i)
!!$             call add_Variable(region,-1     ,'RWET_'//names(itracer)(1:3), 'wet radius of '//names(itracer)(1:3), 'm',3,index_aer,status  )
!!$          end do
!!$          ! dry radius only for soluble, as rwet=rdry for insoluble
!!$          do i=1,nsol
!!$             itracer=mode_start(i)
!!$
!!$             call add_Variable(region,-1     ,'RDRY_'//names(itracer)(1:3), 'dry radius  of '//names(itracer)(1:3), 'm',3,index_aer,status  )
!!$             
!!$          end do

          ! 2d data                    
          ! For now it is exactly the same as in AEROCOM-output, mostly due to optics depending on 
          ! output indices now in place
          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(od550soa     )%mf  = metafields( -1     , 'od550soa     ', 'SOA AOD@550nm                      ' , '1         ', '',     'atmosphere_optical_thickness_due_to_particulate_secondary_organic_aerosol_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         ' )
             mixf(region)%f2d( tr2d_1)%mf  = metafields( iso4nus, 'SO4NUS       ',      'concentration of tracer 01' ,                 'kg m-3', '', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_2)%mf  = metafields( iso4ais, 'SO4AIS       ',      'concentration of tracer 02' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_3)%mf  = metafields( iso4acs, 'SO4ACS       ',      'concentration of tracer 03' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_4)%mf  = metafields( iso4cos, 'SO4COS       ',      'concentration of tracer 04' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_5)%mf  = metafields( ibcais , 'BCAIS        ',      'concentration of tracer 05' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_6)%mf  = metafields( ibcacs , 'BCACS        ',      'concentration of tracer 06' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_7)%mf  = metafields( ibccos , 'BCCOS        ',      'concentration of tracer 07' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_8)%mf  = metafields( ibcaii , 'BCAII        ',      'concentration of tracer 08' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d( tr2d_9)%mf  = metafields( ipomais, 'POMAIS       ',      'concentration of tracer 09' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_10)%mf  = metafields( ipomacs, 'POMACS       ',      'concentration of tracer 10' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_11)%mf  = metafields( ipomcos, 'POMCOS       ',      'concentration of tracer 11' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_12)%mf  = metafields( ipomaii, 'POMAII       ',      'concentration of tracer 12' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_13)%mf  = metafields( issacs , 'SSACS        ',      'concentration of tracer 13' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_14)%mf  = metafields( isscos , 'SSCOS        ',      'concentration of tracer 14' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_15)%mf  = metafields( iduacs , 'DUACS        ',      'concentration of tracer 15' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_16)%mf  = metafields( iducos , 'DUCOS        ',      'concentration of tracer 16' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_17)%mf  = metafields( iduaci , 'DUACI        ',      'concentration of tracer 17' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_18)%mf  = metafields( iducoi , 'DUCOI        ',      'concentration of tracer 18' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_19)%mf  = metafields( isoanus, 'SOANUS       ',      'concentration of tracer 19' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_20)%mf  = metafields( isoaais, 'SOAAIS       ',      'concentration of tracer 20' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_21)%mf  = metafields( isoaacs, 'SOAACS       ',      'concentration of tracer 21' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_22)%mf  = metafields( isoacos, 'SOACOS       ',      'concentration of tracer 22' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )
             mixf(region)%f2d(tr2d_23)%mf  = metafields( isoaaii, 'SOAAII       ',      'concentration of tracer 23' ,                 'kg m-3','', 'concentration_of_tracer_dry_aerosol_in_air'  )

!!$
             mixf(region)%f2d(cc01)%mf  = metafields( inus_n , 'conc_NUS', 'number concentration of mode 01' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(cc02)%mf  = metafields( iais_n , 'conc_AIS', 'number concentration of mode 02' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(cc03)%mf  = metafields( iacs_n , 'conc_ACS', 'number concentration of mode 03' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(cc04)%mf  = metafields( icos_n , 'conc_COS', 'number concentration of mode 04' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(cc05)%mf  = metafields( iaii_n , 'conc_AII', 'number concentration of mode 05' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(cc06)%mf  = metafields( iaci_n , 'conc_ACI', 'number concentration of mode 06' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(cc07)%mf  = metafields( icoi_n , 'conc_COI', 'number concentration of mode 07' ,  'm-3'    ,'', 'number_concentration_of_ambient_aerosol_in_air'  )
             mixf(region)%f2d(h2o1)%mf  = metafields(  -1    , 'aerh2o01 ', 'conc of aerosol water of mode 01' ,  'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air'  )
             mixf(region)%f2d(h2o2)%mf  = metafields(  -1    , 'aerh2o02 ', 'conc of aerosol water of mode 02' ,  'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air'  )
             mixf(region)%f2d(h2o3)%mf  = metafields(  -1    , 'aerh2o03 ', 'conc of aerosol water of mode 03' ,  'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air'  )
             mixf(region)%f2d(h2o4)%mf  = metafields(  -1    , 'aerh2o04 ', 'conc of aerosol water of mode 04' ,  'kg m-3','', 'concentration_of_water_in_ambient_aerosol_in_air'  )
             mixf(region)%f2d(rw01)%mf  = metafields( -1     , 'RWETnus      ', 'rwet nus               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rw02)%mf  = metafields( -1     , 'RWETais      ', 'rwet ais               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rw03)%mf  = metafields( -1     , 'RWETacs      ', 'rwet acs               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rw04)%mf  = metafields( -1     , 'RWETcos      ', 'rwet cos               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rw05)%mf  = metafields( -1     , 'RWETaii      ', 'rwet aii               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rw06)%mf  = metafields( -1     , 'RWETaci      ', 'rwet aci               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rw07)%mf  = metafields( -1     , 'RWETcoi      ', 'rwet coi               ' ,      'm','', 'wet_radius'  )
             mixf(region)%f2d(rd01)%mf  = metafields( -1     , 'RDRYnus      ', 'rdry nus               ' ,      'm','', 'dry radius'  )
             mixf(region)%f2d(rd02)%mf  = metafields( -1     , 'RDRYais      ', 'rdry ais               ' ,      'm','', 'dry radius'  )
             mixf(region)%f2d(rd03)%mf  = metafields( -1     , 'RDRYacs      ', 'rdry acs               ' ,      'm','', 'dry radius'  )
             mixf(region)%f2d(rd04)%mf  = metafields( -1     , 'RDRYcos      ', 'rdry cos               ' ,      'm','', 'dry radius'  )
             mixf(region)%f2d(emiterp)%mf  = metafields( -1     , 'emiterp      ', 'emiterp             ' ,      'kgm-2s-1','', 'emission_of_terpene'  )
             mixf(region)%f2d(emiisop)%mf  = metafields( -1     , 'emiisop      ', 'emiisop              ' ,      'kgm-2s-1','', 'emission_of_isoprene'  )
!!$             mixf(region)%f2d(rd05)%mf  = metafields( -1     , 'RDRYaii      ', 'rdry                   ' ,      'm','', 'precipitation_flux'  )
!!$             mixf(region)%f2d(rd06)%mf  = metafields( -1     , 'RDRYaci      ', 'rdry                   ' ,      'm','', 'precipitation_flux'  )
!!$             mixf(region)%f2d(rd07)%mf  = metafields( -1     , 'RDRYcoi      ', 'rdry                   ' ,      'm','', 'precipitation_flux'  )
!!$             !mixf(region)%f2d(31)%mf  = metafields( ino3_a , 'mmrno3_a    ', 'mmr of nitrate aerosol' ,  'kg m-3', 'mass_fraction_of_nitrate_dry_aerosol_in_air'  )
!!$             !mixf(region)%f2d(32)%mf  = metafields( inh4   , 'mmrnh4      ', 'mmr of ammonium' ,  'kg m-3', 'mass_fraction_of_ammonium_dry_aerosol_in_air'  )
!!$             !mixf(region)%f2d(33)%mf  = metafields( imsa   , 'mmrmsa      ', 'mmr of methane sulfonic acid' ,  'kg m-3', 'mass_fraction_of_methane_sulfonic_acid_in_air'  )
          mixf(region)%f2d(sconcsoa     )%mf  = metafields( -1     , 'sconcsoa     ', 'Surface Concentration SOA          ' , 'kg m-3    ', '',     'mass_concentration_of_secondary_particulate_organic_matter_dry_aerosol_in_air'  )
          mixf(region)%f2d(loadsoa      )%mf  = metafields( -1     , 'loadsoa      ', 'Load of SOA                        ' , 'kg m-2    ', '',     'atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol'  )
          mixf(region)%f2d(drysoa       )%mf  = metafields( -1     , 'drysoa       ', 'Dry Deposition of SOA              ' , 'kg m-2 s-1', 'down', &
                           'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_dry_deposition'  )
          mixf(region)%f2d(wetsoa     )%mf = metafields( -1, 'wetsoa         ', 'Wet Deposition of SOA                ', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_wet_deposition')
          mixf(region)%f2d(p_svoc2d   )%mf = metafields( -1, 'p_svoc2D       ', 'Column integral of SVOC production   ', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_wet_deposition')
          mixf(region)%f2d(p_elvoc2d  )%mf = metafields( -1, 'p_elvoc2D      ', 'Column integral of ELVOC production  ', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_secondary_particulate_organic_matter_dry_aerosol_due_to_wet_deposition')
          mixf(region)%f2d(p_gas_so42d)%mf = metafields( -1, 'prod_gas_so42D ', 'Column integral of gas SO4 production', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_gas_phase_sulfate'  )
          mixf(region)%f2d(p_liq_so42d)%mf = metafields( -1, 'prod_liq_so42D ', 'Column integral of liq SO4 production', 'kg m-2 s-1', 'down', 'tendency_of_atmosphere_mass_content_of_liquid_phase_sulfate'  )

          ! set global dimensions of fields for netcdf definitions
          mixf (region)%nlon      = imr
          mixf (region)%nlat      = jmr
          mixf (region)%nlev      = lmr

          !allocate space for lon, lat, lev information
          allocate(mixf (region)%lon(imr))
          allocate(mixf (region)%lat(jmr))
          allocate(mixf (region)%lev(lmr))

          ! save the lat and lon data for use in output
          mixf (region)%lat=global_lli(region)%lat_deg
          mixf (region)%lon=global_lli(region)%lon_deg
          mixf (region)%lev=(/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34/)

          ! intermediate structures for budgets
          allocate(drydepos(region)%f2dslast(i1:i2,j1:j2,8))
          allocate(wetdepos(region)%f2dslast(i1:i2,j1:j2,8))
          allocate(emission(region)%f2dslast(i1:i2,j1:j2,9))! pom, bc, su, so2, dms, ss, dust, terp, isop 

          ! 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
       ! ----------------
       ! open NetCDF file  (mixf)
       ! ----------------

       call tau2date(itau,idater)

       ! time (in days since 2001-01-01 00:00) 

       !Change to from beginning of a run (idatei)
       !time_reftime6=sdate_simulation
       call date2tau( time_reftime6, itauref )
       
       select case (trim(gen_freq)) 
         case ('hourly')
           write(idates,  '(i4,3i2.2)') idater(1), idater(2), idater(3), idater(4)+1
           !write(idates,  '(i4,2i2.2)') idater(1), idater(2), idater(3)
           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 General output frequency;")'); call goErr
       end select

       call date2tau( idater, itaucur )
       !move the timestamp in middle of the average period
       if (trim(gen_freq)=='hourly') then
          ! if average period is hour move the timestamp 30 min back
          time_shift=30*60
       else if (trim(gen_freq)=='daily') then
          ! if average period is hour move the timestamp 30 min back
          time_shift = 12*60*60
       else if (trim(gen_freq)=='monthly') then
          ! assume 30 day months
          time_shift=15*60*60*24
       end if
       reftime = (itaucur - itauref + time_shift) / 3600.00 !86400.
       !call date2tau(
       !WRITE(str_reftime,'("days since ",i4,"-",i2,"-",i2," ",i2,":",i2)') idatei(1),idatei(2),idatei(3),idatei(4),idatei(5)    
       !WRITE(str_reftime,'("hours since ",i4,"-",i2,"-",i2," ",i2,":",i2)') idatei(1),idatei(2),idatei(3),idatei(4),idatei(5)
       str_reftime=time_units
       ! Changed file name convention to General Control 2015
       mixf(region)%fname = trim(outdir)//'/'//&
            trim(f_dataid)//'_'//&
            trim(f_model) //'_'//&
            trim(gen_exper)//'_'//&
            trim(idates)  //'_'//&
            trim(gen_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_general_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 BACCHUS', status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit, MDF_GLOBAL, 'source',              'TM5-MP revision xxx in SVN repository, including SOA description', 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'     ,        'Tommi Bergman: tommi.bergman@knmi.nl ; 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',          'BACCHUS', 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, 'lon',     imr, lon_dimid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Def_Dim( mixf(region)%funit, 'lat',      jmr, lat_dimid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       !write(3000,*)lmr
       call MDF_Def_Dim( mixf(region)%funit, 'lev',        34, lev_dimid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       !Unlimited time causes a crash in the parallel NETCDF, when unlimited dimension is increased in the file
       !4.4.0 may correct this, but for  now netcdf v4.4.0 on cca is not working
       call MDF_Def_Dim( mixf(region)%funit, 'time',       1, mixf(region)%timeid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! define dimension variables 
       ! longitude
       call MDF_Def_Var( mixf(region)%funit, 'lon', MDF_DOUBLE, &
            (/lon_dimid/), mixf(region)%lonid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit) 
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'units',          'degrees_east', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'axis',          'X', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'long_name',          'longitude', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%lonid , 'standard_name',          'longitude', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! Write out the longitudes
       call MDF_Put_Var( mixf(region)%funit, mixf(region)%lonid, mixf(region)%lon, status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! define latitude variable
       call MDF_Def_Var( mixf(region)%funit, 'lat', MDF_DOUBLE, &
            (/lat_dimid/), mixf(region)%latid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit) 

       !write out the latitude to variable
       call MDF_Put_Var( mixf(region)%funit, mixf(region)%latid, mixf(region)%lat, status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'units',          'degrees_north', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'axis',          'Y', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'long_name',          'latitude', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%latid , 'standard_name',          'latitude', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! lev
       call MDF_Def_Var( mixf(region)%funit, 'lev' , MDF_DOUBLE, &
            (/lev_dimid/), mixf(region)%levid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit) 
       !call MDF_Put_Var( mixf(region)%funit, mixf(region)%levid,  (/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34/), status)
      call MDF_Put_Var( mixf(region)%funit, mixf(region)%levid,  mixf(region)%lev, status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! time
       call MDF_Def_Var( mixf(region)%funit, 'time', MDF_DOUBLE, &
            (/mixf(region)%timeid/), mixf(region)%time_varid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit) 
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%time_varid , 'units',          str_reftime, status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       !call MDF_Put_Att( mixf(region)%funit,mixf(region)%timeid , 'long_name',          'time', status)
       !IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Att( mixf(region)%funit,mixf(region)%time_varid , 'standard_name',          'time', status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)


       ! define variables
       !do i = 1, ntrace!r_3d 
       do i = 1, n_3d_vars!ntracer_3d 
          !write(1234,*)i,trim(mixf(region)%f3d(i)%mf%vname)
          call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f3d(i)%mf%vname), MDF_FLOAT, &
               (/mixf(region)%lonid, mixf(region)%latid, mixf(region)%levid,  mixf(region)%timeid/), 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 
!          write(1234,*)i,trim(mixf(region)%f2d(i)%mf%vname)
          call MDF_Def_Var( mixf(region)%funit, trim(mixf(region)%f2d(i)%mf%vname), MDF_FLOAT, &
               (/mixf(region)%lonid, mixf(region)%latid, mixf(region)%timeid/), 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_Def_Dim( mixf(region)%funit, 'kappa',     1, kappa_dimid, status )
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       call  MDF_Def_Var( mixf(region)%funit, 'year', MDF_FLOAT, &
               (/kappa_dimid/),year_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'month', MDF_FLOAT, &
               (/kappa_dimid/),month_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'kSO4', MDF_FLOAT, &
               (/kappa_dimid/),kso4_varid, status )
       !status = NF90_Def_Var( mixf(region)%funit, 'kSO4', MDF_FLOAT, &
       !     kso4_varid )
       call  MDF_Def_Var( mixf(region)%funit, 'kPOM', MDF_FLOAT, &
               (/kappa_dimid/), kpom_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'kSOA', MDF_FLOAT, &
               (/kappa_dimid/), ksoa_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'kBC', MDF_FLOAT, &
               (/kappa_dimid/), kbc_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'kDU', MDF_FLOAT, &
               (/kappa_dimid/), kdu_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'kSS', MDF_FLOAT, &
               (/kappa_dimid/), kss_varid, status )

       call  MDF_Def_Var( mixf(region)%funit, 'kNO3', MDF_FLOAT, &
               (/kappa_dimid/), kno3_varid, status )
       call  MDF_Def_Var( mixf(region)%funit, 'kNa2SO4', MDF_FLOAT, &
               (/kappa_dimid/), kna2so4_varid, status )

       call  MDF_Def_Var( mixf(region)%funit, 'kMSA', MDF_FLOAT, &
               (/kappa_dimid/), kmsa_varid, status )

       call MDF_EndDef( mixf(region)%funit, status )

       IF_NOTOK_MDF(fid=mixf(region)%funit)

       ! KAPPA values written only once
       ! needed for BACCHUS intercomparison...

       call MDF_Put_Var( mixf(region)%funit, year_varid, (/idater(1)/), status)
       call MDF_Put_Var( mixf(region)%funit, month_varid, (/idater(2)/), status)
       call MDF_Put_Var( mixf(region)%funit, kso4_varid, (/kap_su/), status)
       !status= NF90_Put_Var( mixf(region)%funit, kso4_varid, kap_su)
       !IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, kpom_varid, (/kap_pom/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, kbc_varid, (/kap_bc/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, ksoa_varid, (/kap_soa/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, kss_varid, (/kap_ss/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, kdu_varid, (/kap_du/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, kno3_varid, (/kap_no3/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)

       call MDF_Put_Var( mixf(region)%funit, kna2so4_varid, (/kap_na2so4/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)
       call MDF_Put_Var( mixf(region)%funit, kmsa_varid, (/kap_msa/), status)
       IF_NOTOK_MDF(fid=mixf(region)%funit)





    end do regionloop ! nregions


    call goLabel() ; status = 0

  end  subroutine output_general_init
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_GENERAL_DONE
  !
  ! !DESCRIPTION: Free parameters. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_general_done(status, iregion)
    !
    ! !USES:
    use chem_param,    only : ntracet,ntrace
    !
    ! !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_general_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( gas_output)
       deallocate( mixf (region)%f3d )
       deallocate( mixf (region)%f2d )
       deallocate( drydepos(region)%f2dslast )
       deallocate( wetdepos(region)%f2dslast )
       deallocate( emission(region)%f2dslast )

      

    end do regionloop

    if (nCCNdiag) then
       deallocate(ind_CCN)
    !   deallocate(N_rad)
    !   deallocate(HG_rad)
    end if
    call goLabel() ; status = 0

  end subroutine output_general_done
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  output_general_write
  !
  ! !DESCRIPTION: This routines builds the average by dividing accumulated
  !               data by stack counter. The results are written to file. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_general_write(region, status)
    !
    ! !USES:
    !
    use chem_param,    only : ntracet,ntrace
    use dims,only:itau
    use datetime,      only : tau2date, date2tau

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

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

    ! Time definitions for General

    real                 :: reftime
    integer              :: itauref, time_shift

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


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

    ! -------------
    ! write to file 
    ! -------------
    ! Ncfile needs a timestep
    
    ! define the reference time
    ! time (for now in days since 2001-01-01 00:00) 
    call date2tau( time_reftime6, itauref )

    ! calculate reftime as fractional days from itauref, hence division by 86400 (24h*3600s)
    !call date2tau( idater, itaucur )

    if (trim(gen_freq)=='hourly') then
       ! if average period is hour move the timestamp 30 min forward because current itau is 
       ! at the beginning of the timestep
       time_shift=30.0*60.0
    else if (trim(gen_freq)=='daily') then
       ! if average period is hour move the timestamp 30 min back
       time_shift = 12*60*60
    else if (trim(gen_freq)=='monthly') then
       ! assume 30 day months
       time_shift=15*60*60*24
    end if

    !reftime = (itau - itauref-time_shift) / 3600 !86400. in hours instead of days
    reftime = (itau - itauref+time_shift) / 3600 !86400. in hours instead of days
    !write (1111,*)reftime,itau,itauref,itau-itauref, itau-itauref-time_shift,time_shift
    reftime =reftime-0.5
    !reftime = (itau - itauref) / 86400.    
    !Add time stamp to current
    !Currently only allows 1 step per file, needs to be extended

    do i=1,ntracer_2d
       call MDF_Put_Var( mixf(region)%funit, mixf(region)%f2d(i)%varid, mixf(region)%f2d(i)%field(i1:i2,j1:j2), status, start=(/i1,j1,1/), count=(/imr,jmr,1/) )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
    end do
    do i=1,n_3d_vars!ntracer_3d
       !write(1111,*) mixf(region)%f3d(i)%mf%vname,imr,jmr,lmr, mixf(region)%nlev,mixf(region)%nlat,mixf(region)%nlon
       call MDF_Put_Var( mixf(region)%funit, mixf(region)%f3d(i)%varid, mixf(region)%f3d(i)%field(i1:i2,j1:j2,1:lmr), status, start=(/i1,j1,1,1/), count=(/imr,jmr,lmr,1/) )
       IF_NOTOK_MDF(fid=mixf(region)%funit)
    end do
    call MDF_Var_Par_Access( mixf(region)%funit, mixf(region)%time_varid, MDF_INDEPENDENT, status )
    IF_NOTOK_MDF(fid=mixf(region)%funit)
    
    ! Write current reftime
    call MDF_Put_Var( mixf(region)%funit, mixf(region)%time_varid,(/reftime/), status, start=(/1/),count=(/1/))
    IF_NOTOK_MDF(fid=mixf(region)%funit)
    
    call MDF_Close( mixf(region)%funit, status )
    IF_NOTOK_MDF(fid=mixf(region)%funit)
    
   

    call goLabel() ; status = 0

  end subroutine output_general_write
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    OUTPUT_GENERAL_STEP
  !
  ! !DESCRIPTION: This is the accumulation routine. It is called following 
  !               the user specification general.dhour in rc-file. It 
  !               evaluates the various diagnostics and does summing. 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine output_general_step( region, dhour, status )
    !
    ! !USES:
    !
    use GO         , only : TDate, NewDate, rTotal, operator(-), GO_Timer_Def, GO_Timer_End, GO_Timer_Start
    use Grid       , only : FPressure
    use binas,       only : rgas, rol,pi
    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 chemistry,   only : d_nucle, m_nucle_su, m_nucle_soa, growth1_2, coag_sink_nuc, cond1_soa, cond1_su
    use chemistry,   only : cond2_soa, cond3_soa, cond4_soa, cond5_soa
    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,isoanus
    use chem_param,  only : isoaais, isoaacs, isoacos, isoaaii
    use chem_param,  only : imsa, inh4
    use chem_param,  only : inus_n,iacs_n,icos_n,iais_n,iaii_n,iaci_n,icoi_n
    use chem_param,  only : ntracet,ntrace,xms,xmso4,xmc,xmbc,xmpom
    use chem_param,  only : mode_end_pom,mode_end_so4,mode_end_ss,mode_end_bc,mode_end_dust,mode_end_soa
    use mo_aero_m7,  only : nsol,nmod,dnacl,doc,dh2so4,dbc,ddust
    use optics,      only : optics_aop_get
    use m7_data,     only : h2o_mode, rw_mode, rwd_mode
    use ebischeme,   only : diag_prod
    use mo_aero_m7,  only : cmedr2mmedr,sigma
   !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    integer, intent(in)  :: dhour   ! this is general.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_general_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, iSat
    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
    integer  :: rwetcounter,h2ocounter,rdrycounter
    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

    real, dimension(:,:,:,:),   allocatable :: CCN 
    real, dimension(:,:,:,:),   allocatable :: dry_rad

    real, dimension(:,:,:,:),   allocatable :: Kappa_r 
    real, dimension(:,:,:,:),   allocatable :: N_r

    integer :: ncontr, lvec, lct, l, indoffset, nwl,nn,irw,irwd
    integer :: nadd, nadd_max, indoffset_stat, indabs
    integer :: iadd
    real    :: z,sizelimit,hr2,rad
    real    :: no3fac, wght, dwght, tabs
    real,dimension(7)    :: modfrac,lnsigma
    integer,dimension(6) :: idate_f

    type(TDate)          ::  t, t0
    real                 ::  time 
    real                 ::  perm3_to_percm3,kg_to_ug
    integer              ::  d_ind
    ! --- 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


    ! link the datastructures for transported and non transported variables and 
    ! Gridbox area
    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
   
    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)

    allocate( CCN(nSat,i1:ie,j1:je,lmr))
    allocate( dry_rad(nsol,i1:ie,j1:je,lmr))

    allocate( Kappa_r(3,i1:ie,j1:je,lmr))
    allocate( N_r(3,i1:ie,j1:je,lmr))

    call CCN_Diag(CCN, rm(i1:ie,j1:je,1:lmr,1:ntracet), dry_rad, temper_dat(region)%data(i1:i2,j1:j2,1:lmr))
    call Kappa_diag(CCN, rm(i1:ie,j1:je,1:lmr,1:ntracet), dry_rad, temper_dat(region)%data(i1:i2,j1:j2,1:lmr),N_r,Kappa_r)

    ! 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
   
    ! ----------------------
    ! 3D meteo fields
    ! ----------------------
    
    !temp=1
    !hus=2
    !airmass=3
    !pressure=4


    ! 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
    
    do k=1,lmr
       do j=j1,j2
          do i=i1,i2
             dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * Rgas) 
!!$             mixf(region)%f3d(p_gas_so4)%field(i,j,k)= mixf(region)%f3d(p_gas_so4)%field(i,j,k)+diag_prod(region)%prod(i,j,k,1) / dxyp(j)
!!$             mixf(region)%f3d(p_liq_so4)%field(i,j,k)= mixf(region)%f3d(p_liq_so4)%field(i,j,k)+diag_prod(region)%prod(i,j,k,2) / dxyp(j)
             mixf(region)%f2d(p_gas_so42d)%field(i,j)= mixf(region)%f2d(p_gas_so42d)%field(i,j)+diag_prod(region)%prod(i,j,k,1) / dxyp(j)
             mixf(region)%f2d(p_liq_so42d)%field(i,j)= mixf(region)%f2d(p_liq_so42d)%field(i,j)+diag_prod(region)%prod(i,j,k,2) / dxyp(j)
!!$            
!!$             mixf(region)%f3d(p_elvoc)%field(i,j,k)= mixf(region)%f3d(p_elvoc)%field(i,j,k)+diag_prod(region)%prod(i,j,k,3) / dxyp(j)! divide by 2 to calculate mean over 2 chemistry timesteps in one model timestep
!!$             mixf(region)%f3d(p_svoc)%field(i,j,k)= mixf(region)%f3d(p_svoc)%field(i,j,k)+diag_prod(region)%prod(i,j,k,4) / dxyp(j)
!!$
!!$             mixf(region)%f3d(p_elohterp)%field(i,j,k)= mixf(region)%f3d(p_elohterp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,5) / dxyp(j)
!!$             mixf(region)%f3d(p_elo3terp)%field(i,j,k)= mixf(region)%f3d(p_elo3terp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,6) / dxyp(j)
!!$             mixf(region)%f3d(p_elohisop)%field(i,j,k)= mixf(region)%f3d(p_elohisop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,7) / dxyp(j)
!!$             mixf(region)%f3d(p_elo3isop)%field(i,j,k)= mixf(region)%f3d(p_elo3isop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,8) / dxyp(j)
!!$             mixf(region)%f3d(p_svohterp)%field(i,j,k)= mixf(region)%f3d(p_svohterp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,9) / dxyp(j)
!!$             mixf(region)%f3d(p_svo3terp)%field(i,j,k)= mixf(region)%f3d(p_svo3terp)%field(i,j,k)+diag_prod(region)%prod(i,j,k,10) / dxyp(j)
!!$             mixf(region)%f3d(p_svohisop)%field(i,j,k)= mixf(region)%f3d(p_svohisop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,11) / dxyp(j)
!!$             mixf(region)%f3d(p_svo3isop)%field(i,j,k)= mixf(region)%f3d(p_svo3isop)%field(i,j,k)+diag_prod(region)%prod(i,j,k,12) / dxyp(j)

             mixf(region)%f2d(p_elvoc2d)%field(i,j)= mixf(region)%f2d(p_elvoc2d)%field(i,j)+diag_prod(region)%prod(i,j,k,3) / dxyp(j)
             mixf(region)%f2d(p_svoc2d)%field(i,j)= mixf(region)%f2d(p_svoc2d)%field(i,j)+diag_prod(region)%prod(i,j,k,4) / dxyp(j)
!!$             mixf(region)%f3d(d_nuc)%field(i,j,k)= mixf(region)%f3d(d_nuc)%field(i,j,k)+d_nucle(i,j,k)*dtime !value is saved #cm-3s-1
!!$             mixf(region)%f3d(m_nuc_su)%field(i,j,k)= mixf(region)%f3d(m_nuc_su)%field(i,j,k)+m_nucle_su(i,j,k)*dtime !value is saved #molec cm-3s-1
!!$             mixf(region)%f3d(m_nuc_soa)%field(i,j,k)= mixf(region)%f3d(m_nuc_soa)%field(i,j,k)+m_nucle_soa(i,j,k)*dtime !value is saved #µg m-3s-1
!!$             mixf(region)%f3d(gr1_2)%field(i,j,k)= mixf(region)%f3d(gr1_2)%field(i,j,k)+growth1_2(i,j,k)*dtime !value is saved #cm-3s-1
!!$             mixf(region)%f3d(coag1)%field(i,j,k)= mixf(region)%f3d(coag1)%field(i,j,k)+coag_sink_nuc(i,j,k)*dtime !value is saved #cm-3s-1
!!$             mixf(region)%f3d(co1_soa)%field(i,j,k)= mixf(region)%f3d(co1_soa)%field(i,j,k)+cond1_soa(i,j,k)*dtime !value is saved # unit?
!!$             mixf(region)%f3d(co2_soa)%field(i,j,k)= mixf(region)%f3d(co2_soa)%field(i,j,k)+cond2_soa(i,j,k)*dtime !value is saved # unit?
!!$             mixf(region)%f3d(co3_soa)%field(i,j,k)= mixf(region)%f3d(co3_soa)%field(i,j,k)+cond3_soa(i,j,k)*dtime !value is saved # unit?
!!$             mixf(region)%f3d(co4_soa)%field(i,j,k)= mixf(region)%f3d(co4_soa)%field(i,j,k)+cond4_soa(i,j,k)*dtime !value is saved # unit?
!!$             mixf(region)%f3d(co5_soa)%field(i,j,k)= mixf(region)%f3d(co5_soa)%field(i,j,k)+cond5_soa(i,j,k)*dtime !value is saved # unit?
!!$             mixf(region)%f3d(co1_su)%field(i,j,k)= mixf(region)%f3d(co1_su)%field(i,j,k)+cond1_su(i,j,k)*dtime !value is saved # unit?

             !BACCHUS pm1 vars
             lnsigma=log(sigma)
             hr2= 0.5*sqrt(2.0)
             sizelimit=1.0*5.e-7! pm1(1.0 um in diam) in radius (m)
             Do imode = 1, nmod
                ! Calculate the median radius of the volume weighted distribution.
                rad = rw_mode(region,imode)%d3(i,j,k) * cmedr2mmedr(imode)
                
                !!             if( rad == 0.0 ) then ! changed to a precision depending value
                if( rad < 100*TINY(1.0) ) then
                   modfrac(imode) = 1.0 ! Should not matter, because if the radius is zero, 
                   ! there are no aerosols. But in principle, it would 
                   ! mean that all aerosols are infinitely small, so 
                   ! they would fit in any PM-class.
                else
                   z = ( log(sizelimit) - log(rad) ) / lnsigma(imode)
                   modfrac(imode) = 0.5 + 0.5 * ERF(z*hr2)
                end if
             end do
             perm3_to_percm3=1.0e-6
             kg_to_ug=1.0e9
             mixf(region)%f3d(mmf1)%field(i,j,k) = mixf(region)%f3d(mmf1)%field(i,j,k) + dtime * &
                  modfrac(1)
             mixf(region)%f3d(mmf2)%field(i,j,k) = mixf(region)%f3d(mmf2)%field(i,j,k) + dtime * &
                  modfrac(2)
             mixf(region)%f3d(mmf3)%field(i,j,k) = mixf(region)%f3d(mmf3)%field(i,j,k) + dtime * &
                  modfrac(3)
             mixf(region)%f3d(mmf4)%field(i,j,k) = mixf(region)%f3d(mmf4)%field(i,j,k) + dtime * &
                  modfrac(4)
             mixf(region)%f3d(so4pm1)%field(i,j,k) = mixf(region)%f3d(so4pm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,iso4nus)*modfrac(1)+rm(i,j,k,iso4ais)*modfrac(2)+&
                           rm(i,j,k,iso4acs)*modfrac(3)+rm(i,j,k,iso4cos)*modfrac(4))*(xms/xmso4)*kg_to_ug
             mixf(region)%f3d(bcpm1)%field(i,j,k) = mixf(region)%f3d(bcpm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,ibcais)*modfrac(2)+&
                           rm(i,j,k,ibcacs)*modfrac(3)+rm(i,j,k,ibccos)*modfrac(4)+rm(i,j,k,ibcaii)*modfrac(5))*kg_to_ug*(xmc/xmbc)
             mixf(region)%f3d(orgpm1)%field(i,j,k) = mixf(region)%f3d(orgpm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,isoanus)*modfrac(1)+(rm(i,j,k,isoaais)+rm(i,j,k,ipomais))*modfrac(2)+&
                           (rm(i,j,k,isoaacs)+rm(i,j,k,ipomacs))*modfrac(3)+(rm(i,j,k,isoacos)+rm(i,j,k,ipomcos))*modfrac(4)+&
                           (rm(i,j,k,isoaaii)+rm(i,j,k,ipomaii))*modfrac(5))*kg_to_ug*(xmc/xmpom)
             mixf(region)%f3d(soapm1)%field(i,j,k) = mixf(region)%f3d(soapm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,isoanus)*modfrac(1)+(rm(i,j,k,isoaais))*modfrac(2)+&
                           (rm(i,j,k,isoaacs))*modfrac(3)+(rm(i,j,k,isoacos))*modfrac(4)+&
                           (rm(i,j,k,isoaaii))*modfrac(5))*kg_to_ug*(xmc/xmpom)
             mixf(region)%f3d(sspm1)%field(i,j,k) = mixf(region)%f3d(sspm1)%field(i,j,k) + dtime * &
                           dens * ( rm(i,j,k,issacs)*modfrac(3)+rm(i,j,k,isscos)*modfrac(4))*kg_to_ug
             mixf(region)%f3d(dupm1)%field(i,j,k) = mixf(region)%f3d(dupm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,iduacs)*modfrac(3)+&
                           rm(i,j,k,iducos)*modfrac(4)+rm(i,j,k,iduaci)*modfrac(6)+rm(i,j,k,iducoi)*modfrac(7))*kg_to_ug
            mixf(region)%f3d(nh4pm1)%field(i,j,k) = mixf(region)%f3d(nh4pm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,inh4)*modfrac(3))*kg_to_ug
                           
            mixf(region)%f3d(no3pm1)%field(i,j,k) = mixf(region)%f3d(no3pm1)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,ino3_a)*modfrac(3))*kg_to_ug
                           
             mixf(region)%f3d(so4)%field(i,j,k) = mixf(region)%f3d(so4)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,iso4nus)+rm(i,j,k,iso4ais)+&
                           rm(i,j,k,iso4acs)+rm(i,j,k,iso4cos))*(xms/xmso4)*kg_to_ug
             mixf(region)%f3d(bc)%field(i,j,k) = mixf(region)%f3d(bc)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,ibcais)+&
                           rm(i,j,k,ibcacs)+rm(i,j,k,ibccos)+rm(i,j,k,ibcaii))*(xmc/xmbc)*kg_to_ug
             mixf(region)%f3d(org)%field(i,j,k) = mixf(region)%f3d(org)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,isoanus)+(rm(i,j,k,isoaais)+rm(i,j,k,ipomais))+&
                           (rm(i,j,k,isoaacs)+rm(i,j,k,ipomacs))+(rm(i,j,k,isoacos)+rm(i,j,k,ipomcos))+&
                           (rm(i,j,k,isoaaii)+rm(i,j,k,ipomaii)))*(xmc/xmpom)*kg_to_ug
             mixf(region)%f3d(ss)%field(i,j,k) = mixf(region)%f3d(ss)%field(i,j,k) + dtime * &
                           dens * ( rm(i,j,k,issacs)+rm(i,j,k,isscos))*kg_to_ug
             mixf(region)%f3d(du)%field(i,j,k) = mixf(region)%f3d(du)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,iduacs)+&
                           rm(i,j,k,iducos)+rm(i,j,k,iduaci)+rm(i,j,k,iducoi))*kg_to_ug
            mixf(region)%f3d(nh4)%field(i,j,k) = mixf(region)%f3d(nh4)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,inh4))*kg_to_ug
                           
            mixf(region)%f3d(no3)%field(i,j,k) = mixf(region)%f3d(no3)%field(i,j,k) + dtime * &
                           dens * (rm(i,j,k,ino3_a))*kg_to_ug
                           

            
            sizelimit=3*5.e-9! pm1(1.0 um in diam) in radius (m)
            Do imode = 1, nmod
                ! Calculate the median radius of the volume weighted distribution.
                rad = rw_mode(region,imode)%d3(i,j,k) 
                if( rad < 100*TINY(1.0) ) then
                   modfrac(imode) = 1.0 ! Should not matter, because if the radius is zero, 
                   ! there are no aerosols. But in principle, it would 
                   ! mean that all aerosols are infinitely small, so 
                   ! they would fit in any PM-class.
                else
                   z = ( log(sizelimit) - log(rad) ) / lnsigma(imode)
                   
                   ! now calculating larger than sizelimit
                   modfrac(imode) = 1-(0.5 + 0.5 * ERF(z*hr2))
                end if
             end do
             mixf(region)%f3d(n3)%field(i,j,k) = mixf(region)%f3d(n3)%field(i,j,k) + dtime * &
                     dens *( rm(i,j,k,inus_n)*modfrac(1)+ rm(i,j,k,iais_n)*modfrac(2)+&
                     rm(i,j,k,iacs_n)*modfrac(3)+ rm(i,j,k,icos_n)*modfrac(4)+&
                     rm(i,j,k,iaii_n)*modfrac(5)+ rm(i,j,k,iaci_n)*modfrac(6)+&
                     rm(i,j,k,icoi_n)*modfrac(7))*perm3_to_percm3
!!$             !50nm
!!$             do d_ind=1,3
!!$                sizelimit=N_diam(d_ind)*5.e-9! pm1(1.0 um in diam) in radius (m)
!!$                Do imode = 1, nmod
!!$                   ! Calculate the median radius of the volume weighted distribution.
!!$                   rad = rw_mode(region,imode)%d3(i,j,k) 
!!$                   if( rad < 100*TINY(1.0) ) then
!!$                      modfrac(imode) = 1.0 ! Should not matter, because if the radius is zero, 
!!$                      ! there are no aerosols. But in principle, it would 
!!$                      ! mean that all aerosols are infinitely small, so 
!!$                      ! they would fit in any PM-class.
!!$                   else
!!$                      z = ( log(sizelimit) - log(rad) ) / lnsigma(imode)
!!$                      
!!$                      ! now calculating larger than sizelimit
!!$                      modfrac(imode) = 1-(0.5 + 0.5 * ERF(z*hr2))
!!$                   end if
!!$                end do
!!$                mixf(region)%f3d(ind_N(d_ind))%field(i,j,k) = mixf(region)%f3d(n3)%field(i,j,k) + dtime * &
!!$                     dens *( rm(i,j,k,inus_n)*modfrac(1)+ rm(i,j,k,iais_n)*modfrac(2)+&
!!$                     rm(i,j,k,iacs_n)*modfrac(3)+ rm(i,j,k,icos_n)*modfrac(4)+&
!!$                     rm(i,j,k,iaii_n)*modfrac(5)+ rm(i,j,k,iaci_n)*modfrac(6)+&
!!$                     rm(i,j,k,icoi_n)*modfrac(7))*perm3_to_percm3
!!$             end do
 
             
             
             ! CCN under given supersaturation
             IF (nCCNdiag) THEN
                DO iSat=1,nSat
                   mixf(region)%f3d(ind_CCN(iSat))%field(i,j,k)= mixf(region)%f3d(ind_CCN(iSat))%field(i,j,k) + dtime * &
                       dens * 1.e-06 * CCN(iSat,i,j,k) ! cm-3
                ENDDO
 
             DO iSat=1,3
                mixf(region)%f3d(ind_HG(iSat))%field(i,j,k)= mixf(region)%f3d(ind_HG(iSat))%field(i,j,k) + dtime * &
                      Kappa_r(iSat,i,j,k) ! cm-3
                mixf(region)%f3d(ind_N(iSat))%field(i,j,k)= mixf(region)%f3d(ind_N(iSat))%field(i,j,k) + dtime * &
                     dens * 1.e-06 * N_r(iSat,i,j,k) ! cm-3
             ENDDO
            ENDIF
             mixf(region)%f3d(ccn02)%field(i,j,k) = mixf(region)%f3d(ccn02)%field(i,j,k) + dtime * &
                       dens * 1.e-06 * CCN(5,i,j,k) ! cm-3 CCN0.2
             mixf(region)%f3d(gph)%field(i,j,k) = mixf(region)%f3d(gph)%field(i,j,k) + gph_dat(region)%data(i,j,k)*dtime 


             ! rwet, rdry, h2o data holder indices go from 1-4/7 
             ! so use counters, could also add modeindex in mixf
             rwetcounter=1
             rdrycounter=1
             h2ocounter=1
             
             do nn=1,n_3d_vars
                
                if (mixf(region)%f3d(nn)%mf%itm5>0)then
                   ! transported variables in rm
                   if (mixf(region)%f3d(nn)%mf%itm5<=ntracet) then
                      mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * &
                           dens * rm(i,j,k,mixf(region)%f3d(nn)%mf%itm5)
                   else !non transported variables in rmc
                       mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * &
                           dens * rmc(i,j,k,mixf(region)%f3d(nn)%mf%itm5)
                    end if
                    
                 else  ! not a tracer variable
                    ! counter to go through 7 modes when outputing rwet
                    ! requires that rwet variables are in right order
                    ! Wet radius of particles
                   if ( mixf(region)%f3d(nn)%mf%vname(1:4)=='RWET')then

                      mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * &

                           rw_mode(region,rwetcounter)%d3(i,j,k)
                      rwetcounter=rwetcounter+1
                   end if
                   ! aerosol water
                   if ( mixf(region)%f3d(nn)%mf%vname(1:6)=='aerh2o')then
                      mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * &
                           dens*h2o_mode(region,h2ocounter)%d3(i,j,k)
                      
                      h2ocounter=h2ocounter+1
                   end if
  
                   !dry radius of soluble modes
                   if ( mixf(region)%f3d(nn)%mf%vname(1:4)=='RDRY')then
                      mixf(region)%f3d(nn)%field(i,j,k) = mixf(region)%f3d(nn)%field(i,j,k) + dtime * &
                           rwd_mode(region,rdrycounter)%d3(i,j,k)
                      rdrycounter=rdrycounter+1
                   end if

                end if

             end do
          end do
       end do
    end do

    ! set fields summed in chemistry to 0 after writing
    d_nucle = 0.
    m_nucle_su = 0.
    m_nucle_soa = 0.
    diag_prod(region)%prod(i1:i2,j1:j2,:,1)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,2)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,3)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,4)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,5)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,6)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,7)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,8)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,9)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,10)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,11)=0.0
    diag_prod(region)%prod(i1:i2,j1:j2,:,12)=0.0
    coag_sink_nuc = 0.
    growth1_2 = 0.
    cond1_soa = 0.
    cond2_soa = 0.
    cond3_soa = 0.
    cond4_soa = 0.
    cond5_soa = 0.
    cond1_su = 0.

    ! ----------------------
    ! 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]
          ! ----------------------
          !TB added nucleation mode oc
          ! 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)) 
          ! SOA: (NUS + AIS + ACS + COS + AII)
          mixf(region)%f2d(sconcsoa)%field(i,j) = mixf(region)%f2d(sconcsoa)%field(i,j) + dtime * & 
               dens * (rm(i,j,k,iSOAnus) + rm(i,j,k,iSOAais) + rm(i,j,k,iSOAacs) + rm(i,j,k,iSOAcos) + rm(i,j,k,iSOAaii))
          ! 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)

          ! Mass concentrations
          do nn=68,90
             mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) + dtime * &
                  dens * rm(i,j,k,mixf(region)%f2d(nn)%mf%itm5)
          end do
          !number concentrations
          do nn=91,97
             mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) + dtime * &
                  dens * rm(i,j,k,mixf(region)%f2d(nn)%mf%itm5)
          end do
          !h2o
          h2ocounter=0
          do nn=98,101
             h2ocounter=h2ocounter+1
             mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) + dtime * &
                  dens *  h2o_mode(region,h2ocounter)%d3(i,j,k)!rm(i,j,k,mixf(region)%f2d(nn)%mf%itm5)
          end do

          ! RWET for all modes  (RDRY=RWET for non-soluble)
          ! 
          irw=0
          do nn=102,108
          !
             irw=irw+1
             mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) +dtime *rw_mode(region,irw)%d3(i,j,k)
          end do
          ! RDRY for soluble modes
          irwd=0
          do nn=109,112
             irwd=irwd+1
             mixf(region)%f2d(nn)%field(i,j) = mixf(region)%f2d(nn)%field(i,j) +dtime *rwd_mode(region,irwd)%d3(i,j,k)
          end do



          ! ----------------------
          ! Load in [kg m-2]
          ! ----------------------
          do k = 1, lmr 
             !do ll=1,len(load)

             ! POM: (AIS + ACS + COS + AII)
             load_tmp=0.0
             do nn=1,7
             
                if (mode_end_pom(nn)>0) then
                   load_tmp = load_tmp+ rm(i,j,k,mode_end_pom(nn))
                   !(rm(i,j,k,iPOMnus) + rm(i,j,k,iPOMais) + rm(i,j,k,iPOMacs) + rm(i,j,k,iPOMcos) + rm(i,j,k,iPOMaii))
                end if
             end do
             mixf(region)%f2d(loadoa)%field(i,j) = mixf(region)%f2d(loadoa)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! SOA: (NUS + AIS + ACS + COS + AII)
             load_tmp=0.0
             do nn=1,7
             
                if (mode_end_soa(nn)>0) then
                   load_tmp = load_tmp+ rm(i,j,k,mode_end_soa(nn))
                ENDIF
             end do
             mixf(region)%f2d(loadsoa)%field(i,j) = mixf(region)%f2d(loadsoa)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! BC: (AIS + ACS + COS + AII)
             load_tmp=0.0
             do nn=1,7
             
                if (mode_end_bc(nn)>0) then
                   load_tmp = load_tmp+ rm(i,j,k,mode_end_bc(nn))
                   !load_tmp = (rm(i,j,k,iBCais) + rm(i,j,k,iBCacs) + rm(i,j,k,iBCcos) + rm(i,j,k,iBCaii))
                end if
             end do
             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=0.0
             do nn=1,7
             
                if (mode_end_so4(nn)>0) then
                    load_tmp = load_tmp+ rm(i,j,k,mode_end_so4(nn))
                    !load_tmp = (rm(i,j,k,iSO4nus) + rm(i,j,k,iSO4ais) + rm(i,j,k,iSO4acs) + rm(i,j,k,iSO4cos))
                end if
             end do
             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=0.0
             do nn=1,7
             
                if (mode_end_dust(nn)>0) then
                   load_tmp = load_tmp+ rm(i,j,k,mode_end_dust(nn))
                   !load_tmp = (rm(i,j,k,iDUacs) + rm(i,j,k,iDUcos) + rm(i,j,k,iDUaci) + rm(i,j,k,iDUcoi))
                end if
             end do

             mixf(region)%f2d(loaddust)%field(i,j) = mixf(region)%f2d(loaddust)%field(i,j) + load_tmp * dtime / dxyp(j)

             ! Seasalt: (ACS + COS)
             load_tmp=0.0
             do nn=1,7
             
                if (mode_end_ss(nn)>0) then
                   load_tmp = load_tmp+ rm(i,j,k,mode_end_ss(nn))
                   !load_tmp = (rm(i,j,k,iSSacs) + rm(i,j,k,iSScos))
                end if
             end do
             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


!!$

       
 
    call GO_Timer_Start( itim_opt_out, status )
    IF_NOTOK_RETURN(status=1)
    ! ---------------------
    ! 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 GENERAL 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
    ncontr = 12

    ! 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 
    allocate( opt_add( i1:i2, j1:j2, lmr, size(wdep_out), nadd   ) ) ; opt_add = 0.0


    ! ---------------------------------
    ! 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_general: 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_general: wrong wavelength given for general 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



    end do

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

    deallocate( Kappa_r)
    deallocate( N_r)


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

    call goLabel() ; status = 0

  end subroutine output_general_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 : isoacos, isoaaii, isoaacs, isoaais, isoanus
    use chem_param,  only : idms, iso2, ntracet, ntrace,iterp,iisop
    use chem_param,  only : xmso2, xmso4, xmdms, xmpom, xmbc, xmdust, xmnacl,xmterp,xmisop
#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

    ! terpene:
    tmpbud2d = sum(collect4d(:,:,:,iterp ),3) * xmterp
    mixf(region)%f2d(emiterp)%field = tmpbud2d - emission(region)%f2dslast(:,:,8)
    emission(region)%f2dslast(:,:,8) = tmpbud2d

    ! isoprene:
    tmpbud2d = sum(collect4d(:,:,:,iisop ),3) * xmisop
    mixf(region)%f2d(emiisop)%field = tmpbud2d - emission(region)%f2dslast(:,:,9)
    emission(region)%f2dslast(:,:,9) = 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

    ! SOA (NUS + AIS + ACS + COS + AII)
    tmpbud2d = ( collect3d(:,:,iSOAais ) + collect3d(:,:,iSOAacs ) + &
                 collect3d(:,:,iSOAcos ) + collect3d(:,:,iSOAaii ) + &
                 collect3d(:,:,iSOAnus )) * xmpom
    mixf(region)%f2d(drysoa)%field = tmpbud2d - drydepos(region)%f2dslast(:,:,8)
    drydepos(region)%f2dslast(:,:,8) = 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

    ! SOA (NUS + AIS + ACS + COS + AII)
    tmpbud2d = ( sum(collect4d(:,:,:,iSOAais),3) + sum(collect4d(:,:,:,iSOAacs),3) + &
                 sum(collect4d(:,:,:,iSOAcos),3) + sum(collect4d(:,:,:,iSOAaii),3) + &
                 sum(collect4d(:,:,:,iSOAnus),3)) * xmpom
    mixf(region)%f2d(wetsoa)%field = tmpbud2d - wetdepos(region)%f2dslast(:,:,8)
    wetdepos(region)%f2dslast(:,:,8) = 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
subroutine add_variable(region,itm5,varname,longname,unit,dims,n_out,status)
  use chem_param, only: mode_end_so4,mode_end_pom,mode_end_bc,mode_end_ss,mode_end_dust
  
  implicit none
  integer:: itm5
  integer:: fileunit
  integer:: varid
  character(len=*)::varname
  character(len=*)::longname
  character(len=*)::unit
  integer:: dims
  integer:: region
  integer,intent(out)::n_out
  integer,intent(out)::status
  logical :: writeout=.true.
  character(len=*), parameter ::  rname = mname//'/output_aerchemmip_add_variable'
  integer ::i1,i2,j1,j2,jmr,imr,lmr
  
  call goLabel(rname)
  
  !call GO_Timer_Start( itim_addvar, status )
  !IF_NOTOK_RETURN(status=1)
  
  
  !write(gol,'("add_variable1")')
  !call goErr
  
  if( (n_2d_vars+1>ntracer_2d) .or.(n_3d_vars+1>ntracer_3d) ) then
     status=1
     IF_ERROR_RETURN(status=1)
  end if
  call Get_DistGrid( dgrid(1), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
  
  ! define sizes for arrays
  imr=i2-i1+1
  jmr=j2-j1+1
  lmr = levi%nlev
  
  ! 2D variables
  if (dims==2) then
     !init variable
     n_2d_vars=n_2d_vars+1
     mixf(region)%f2d(n_2d_vars)%mf  = metafields( itm5 , varname,longname,unit,'','')
     mixf(region)%f2d(n_2d_vars)%writeout=writeout

     !Later make the allocation here
     !allocate( mixf(region)%f2d(n_2d_vars)%field(i1:i2,j1:j2) )
     mixf(region)%f2d(n_2d_vars)%field(i1:i2,j1:j2)=0.0
     n_out=n_2d_vars
  else if (dims==3) then
     !init variable
     n_3d_vars=n_3d_vars+1
     mixf(region)%f3d(n_3d_vars)%mf  = metafields( itm5 , varname,longname,unit,'','')
     mixf(region)%f3d(n_3d_vars)%writeout=writeout

     !Later make the allocation here
     !allocate( mixf(region)%f3d(n_3d_vars)%field(i1:i2,j1:j2,lmr) )
     mixf(region)%f3d(n_3d_vars)%field(i1:i2,j1:j2,lmr)=0.0
     n_out=n_3d_vars
  end if

  call goLabel()
  status = 0
  !call GO_Timer_End( itim_addvar, status )
  !IF_NOTOK_RETURN(status=1)
end subroutine add_variable

subroutine CCN_Diag(CCN,mass_num,cmr,temper)
  use chem_param,    only : iso4ais, iso4acs, iso4cos, ibcais, ibcacs, ibccos, ipomais, ipomacs
  use chem_param,    only : ipomcos, issacs, isscos, iduacs, iducos, isoaais, isoaacs, isoacos
  use chem_param,    only : iais_n, iacs_n, icos_n, sigma_lognormal
  use chem_param,    only : ino3_a, imsa, nh4no3_factor, nh4no3_density, msa_density
  use chem_param,    only : SurfTen, mode_tracers
  use chem_param,    only : Kap_su, Kap_pom, Kap_soa, Kap_bc, Kap_du, Kap_ss, Kap_na2so4, Kap_no3, Kap_msa
  use chem_param,    only : so4_density, pom_density, soa_density, carbon_density, ss_density, dust_density
  use mo_aero_m7,    only : cmr2ram, nsol, ram2cmr
  use mo_aero_m7,    only : dNaCl, dH2SO4, dNa2SO4, wNaCl, wH2SO4, wNa2SO4, wSO4
  use binas,         only : pi, Rgas, xm_h2o, rol

  implicit none

  integer :: iSat, iMode
  real :: CCN(:,:,:,:), mass_num(:,:,:,:), cmr(:,:,:,:), temper(:,:,:)
  real, allocatable :: Vol(:,:,:), Kappa(:,:,:), A_Koehler(:,:,:)
  real, allocatable :: N_act(:,:,:), rc(:,:,:)
  real, allocatable :: nNa(:,:,:), nCl(:,:,:), nNaCl(:,:,:), nH2SO4(:,:,:), nNa2SO4(:,:,:), nSO4(:,:,:)

  CCN(:,:,:,:) = 0.e0
  cmr(:,:,:,:) = 0.e0

  allocate(Vol(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Kappa(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(A_Koehler(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(N_act(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(rc(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nNa(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nCl(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nSO4(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nNaCl(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nNa2SO4(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nH2SO4(size(CCN,2),size(CCN,3),size(CCN,4)))

  do iMode = 2,nsol
     ! Kappa
     SELECT CASE (iMode)
        CASE(2) ! no sea salt and dust in soluble aitken mode
           nSO4= mass_num(:,:,:,iso4ais)/(1.e-3*wSO4)
           nH2SO4 = nSO4 
           Vol  = ( nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3)*(wH2SO4/wSO4) & !mass_num(:,:,:,iso4ais)/so4_density &
                  + mass_num(:,:,:,ipomais)/pom_density &
                  + mass_num(:,:,:,isoaais)/soa_density &
                  + mass_num(:,:,:,ibcais)/carbon_density) ! m3/gridbox
           Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3)*(wH2SO4/wSO4) &!mass_num(:,:,:,iso4ais)/so4_density &
                  + Kap_pom*mass_num(:,:,:,ipomais)/pom_density &
                  + Kap_soa*mass_num(:,:,:,isoaais)/soa_density &
                  + Kap_bc*mass_num(:,:,:,ibcais)/carbon_density) / Vol
           cmr(iMode,:,:,:) = ram2cmr(iMode)*(Vol/mass_num(:,:,:,iais_n)*0.75/pi)**(1./3.) ! m

        CASE(3)
           nNa = mass_num(:,:,:,issacs)/(1.e-3*wNaCl) ! mol/gridbox
           nCl = nNa
           nSO4= mass_num(:,:,:,iso4acs)/(1.e-3*wSO4)
           nNa2SO4 = MIN(nNa/2.e0, nSO4)
           nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa
           nNaCl = MIN(nCl, nNa)
           nCl = nNaCl ! overshoot nCl is supposed to evaporate
           nH2SO4 = nSO4 - nNa2SO4
           
           Vol  = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                  + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                  + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &
                  + mass_num(:,:,:,ipomacs)/pom_density &
                  + mass_num(:,:,:,isoaacs)/soa_density &
                  + mass_num(:,:,:,ibcacs)/carbon_density &
                  + mass_num(:,:,:,iduacs)/dust_density &
                  + mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density &
                  + mass_num(:,:,:,imsa)/msa_density) ! m3/gridbox
           Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &
                  + Kap_na2so4*nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                  + Kap_pom*mass_num(:,:,:,ipomacs)/pom_density &
                  + Kap_soa*mass_num(:,:,:,isoaacs)/soa_density &
                  + Kap_bc*mass_num(:,:,:,ibcacs)/carbon_density &
                  + Kap_ss*nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                  + Kap_du*mass_num(:,:,:,iduacs)/dust_density &
                  + Kap_no3*mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density &
                  + Kap_msa*mass_num(:,:,:,imsa)/msa_density) / Vol
           cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,iacs_n)*0.75/pi)**(1./3.)*ram2cmr(iMode) ! m

        CASE(4)
           nNa = mass_num(:,:,:,isscos)/(1.e-3*wNaCl) ! mol/gridbox
           nCl = nNa
           nSO4= mass_num(:,:,:,iso4cos)/(1.e-3*wSO4)
           nNa2SO4 = MIN(nNa/2.e0, nSO4)
           nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa
           nNaCl = MIN(nCl, nNa)
           nCl = nNaCl ! overshoot nCl is supposed to evaporate
           nH2SO4 = nSO4 - nNa2SO4
           Vol  = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                  + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                  + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &
                  + mass_num(:,:,:,ipomcos)/pom_density &
                  + mass_num(:,:,:,isoacos)/soa_density &
                  + mass_num(:,:,:,ibccos)/carbon_density &
                  + mass_num(:,:,:,iducos)/dust_density) ! m3/gridbox
           Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &                 
                  + Kap_na2so4*nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                  + Kap_pom*mass_num(:,:,:,ipomcos)/pom_density &
                  + Kap_soa*mass_num(:,:,:,isoacos)/soa_density &
                  + Kap_bc*mass_num(:,:,:,ibccos)/carbon_density &
                  + Kap_ss*nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                  + Kap_du*mass_num(:,:,:,iducos)/dust_density) / Vol
           cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,icos_n)*0.75/pi)**(1./3.)*ram2cmr(iMode) ! m
     END SELECT

     where(Kappa < 0.04) Kappa = 0.04

     A_Koehler = 2.0 * xm_h2o * SurfTen / Rgas / rol / temper

     do iSat=1,nSat
        ! critical radius per mode
        rc = A_Koehler/3.0 * (2.0/SuperSat(iSat)/sqrt(Kappa))**(2.0/3.0)

        ! number of activated particles
        N_act = mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc/cmr(iMode,:,:,:)) / &
                sqrt(2.e0)/log(sigma_lognormal(iMode)) )

        CCN(iSat,:,:,:) = CCN(iSat,:,:,:) + N_act ! #/gridbox
     enddo
  enddo

  deallocate(Vol)
  deallocate(Kappa)
  deallocate(A_Koehler)
  deallocate(rc)
  deallocate(N_act)
end subroutine CCN_Diag

subroutine Kappa_diag(CCN,mass_num,cmr,temper,N_r,Kappa_r)
  use chem_param,    only : iso4ais, iso4acs, iso4cos, ibcais, ibcacs, ibccos, ipomais, ipomacs
  use chem_param,    only : ipomcos, issacs, isscos, iduacs, iducos, isoaais, isoaacs, isoacos
  use chem_param,    only : iais_n, iacs_n, icos_n, sigma_lognormal
  use chem_param,    only : ino3_a, imsa, nh4no3_factor, nh4no3_density, msa_density
  use chem_param,    only : SurfTen, mode_tracers
  use chem_param,    only : Kap_su, Kap_pom, Kap_soa, Kap_bc, Kap_du, Kap_ss, Kap_na2so4, Kap_no3, Kap_msa
  use chem_param,    only : so4_density, pom_density, soa_density, carbon_density, ss_density, dust_density
  use mo_aero_m7,    only : cmr2ram, nsol, ram2cmr
  use mo_aero_m7,    only : dNaCl, dH2SO4, dNa2SO4, wNaCl, wH2SO4, wNa2SO4, wSO4
  use binas,         only : pi, Rgas, xm_h2o, rol

  implicit none

  integer :: iSat, iMode,i_rad
  real :: CCN(:,:,:,:), mass_num(:,:,:,:), cmr(:,:,:,:), temper(:,:,:)
  real :: N_r(:,:,:,:), Kappa_r(:,:,:,:)
  real, allocatable :: Vol(:,:,:), Kappa(:,:,:), A_Koehler(:,:,:)
  real, allocatable :: N_act(:,:,:), rc(:)
  real, allocatable :: nNa(:,:,:), nCl(:,:,:), nNaCl(:,:,:), nH2SO4(:,:,:), nNa2SO4(:,:,:), nSO4(:,:,:)
  real, allocatable :: Volss(:,:,:)
  real, allocatable :: Voldu(:,:,:)
  real, allocatable :: Volsu(:,:,:)
  real, allocatable :: Volbc(:,:,:)
  real, allocatable :: Volpom(:,:,:)
  real, allocatable :: Volsoa(:,:,:)
  real, allocatable :: Volna2so4(:,:,:)
  real, allocatable :: Volmsa(:,:,:)
  real, allocatable :: Volno3(:,:,:)
  real, allocatable :: Volsum(:,:,:)
  real, allocatable :: over_rad_frac(:,:,:)

 ! CCN(:,:,:,:) = 0.e0
  cmr(:,:,:,:) = 0.e0
  N_r(:,:,:,:) = 0.e0
  Kappa_r(:,:,:,:) = 0.e0
  allocate(Vol(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volsu(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volss(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Voldu(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volbc(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volsoa(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volpom(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volna2so4(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volmsa(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volno3(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(Volsum(size(CCN,2),size(CCN,3),size(CCN,4)))

  allocate(over_rad_frac(size(CCN,2),size(CCN,3),size(CCN,4)))

  volsum(:,:,:)=0.e0
  volss(:,:,:)=0.e0
  voldu(:,:,:)=0.e0
  volbc(:,:,:)=0.e0
  volsu(:,:,:)=0.e0
  volpom(:,:,:)=0.e0
  volsoa(:,:,:)=0.e0
  volna2so4(:,:,:)=0.e0
  volmsa(:,:,:)=0.e0
  volno3(:,:,:)=0.e0
  over_rad_frac(:,:,:)=0.e0

!  allocate(N_r(3,size(CCN,2),size(CCN,3),size(CCN,4)))
!  allocate(Kappa_r(size(CCN,2),size(CCN,3),size(CCN,4)))


  allocate(Kappa(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(A_Koehler(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(N_act(size(CCN,2),size(CCN,3),size(CCN,4)))
!  allocate(rc(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(rc(3))
  allocate(nNa(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nCl(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nSO4(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nNaCl(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nNa2SO4(size(CCN,2),size(CCN,3),size(CCN,4)))
  allocate(nH2SO4(size(CCN,2),size(CCN,3),size(CCN,4)))

  !Diam for N
  !BACCHUS defs D=50,80,120nm->r=25,40,60

  rc(1)= 25.0e-9
  rc(2)= 40.0e-9
  rc(3)= 60.0e-9

  do i_rad =1,3

     !neglect nucleation mode
     do iMode = 2,nsol
        ! Kappa
        SELECT CASE (iMode)
        CASE(2) ! no sea salt and dust in soluble aitken mode
           nSO4= mass_num(:,:,:,iso4ais)/(1.e-3*wSO4)
           nH2SO4 = nSO4 

           Vol  = ( nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &! mass_num(:,:,:,iso4ais)/so4_density &
                + mass_num(:,:,:,ipomais)/pom_density &
                + mass_num(:,:,:,isoaais)/soa_density &
                + mass_num(:,:,:,ibcais)/carbon_density) ! m3/gridbox
           Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &!mass_num(:,:,:,iso4ais)/so4_density &
                + Kap_pom*mass_num(:,:,:,ipomais)/pom_density &
                + Kap_soa*mass_num(:,:,:,isoaais)/soa_density &
                + Kap_bc*mass_num(:,:,:,ibcais)/carbon_density) / Vol

           cmr(iMode,:,:,:) = ram2cmr(iMode)*(Vol/mass_num(:,:,:,iais_n)*0.75/pi)**(1./3.) ! m

           !fraction over threshold
           over_rad_frac= 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / &
                sqrt(2.e0)/log(sigma_lognormal(iMode)) )

           volsu= volsu + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) * over_rad_frac           
           volpom = volpom + mass_num(:,:,:,ipomais)/pom_density * over_rad_frac
           volsoa = volsoa + mass_num(:,:,:,isoaais)/soa_density * over_rad_frac
           volbc = volbc + mass_num(:,:,:,ibcais)/carbon_density * over_rad_frac

           N_r(i_rad,:,:,:) = N_r(i_rad,:,:,:)+mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / &
                sqrt(2.e0)/log(sigma_lognormal(iMode)) )
           

        CASE(3)

           nNa = mass_num(:,:,:,issacs)/(1.e-3*wNaCl) ! mol/gridbox
           nCl = nNa
           nSO4= mass_num(:,:,:,iso4acs)/(1.e-3*wSO4)
           nNa2SO4 = MIN(nNa/2.e0, nSO4)
           nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa
           nNaCl = MIN(nCl, nNa)
           nCl = nNaCl ! overshoot nCl is supposed to evaporate
           nH2SO4 = nSO4 - nNa2SO4
           Vol  = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &
                + mass_num(:,:,:,ipomacs)/pom_density &
                + mass_num(:,:,:,isoaacs)/soa_density &
                + mass_num(:,:,:,ibcacs)/carbon_density &
                + mass_num(:,:,:,iduacs)/dust_density &
                + mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density &
                + mass_num(:,:,:,imsa)/msa_density) ! m3/gridbox
           Kappa= ( Kap_su*nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &
                + Kap_na2so4*nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                + Kap_pom*mass_num(:,:,:,ipomacs)/pom_density &
                + Kap_soa*mass_num(:,:,:,isoaacs)/soa_density &
                + Kap_bc*mass_num(:,:,:,ibcacs)/carbon_density &
                + Kap_ss*nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                + Kap_du*mass_num(:,:,:,iduacs)/dust_density &
                + Kap_no3*mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density &
                + Kap_msa*mass_num(:,:,:,imsa)/msa_density) / Vol

           cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,iacs_n)*0.75/pi)**(1./3.)*ram2cmr(iMode) ! m
           over_rad_frac= 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / &
                sqrt(2.e0)/log(sigma_lognormal(iMode)) )
           
           volsu= volsu + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) * over_rad_frac
           volna2so4= volna2so4 + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) * over_rad_frac
            volpom = volpom + mass_num(:,:,:,ipomacs)/pom_density * over_rad_frac
           volsoa = volsoa + mass_num(:,:,:,isoaacs)/soa_density * over_rad_frac
           volbc = volbc + mass_num(:,:,:,ibcacs)/carbon_density * over_rad_frac
           volss =volss + nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) * over_rad_frac
           voldu = voldu + mass_num(:,:,:,iduacs)/dust_density * over_rad_frac
           volno3 = volno3 + mass_num(:,:,:,ino3_a)*nh4no3_factor/nh4no3_density * over_rad_frac
           volmsa = volmsa + mass_num(:,:,:,imsa)/msa_density * over_rad_frac

           N_r(i_rad,:,:,:) = N_r(i_rad,:,:,:)+mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / &
                sqrt(2.e0)/log(sigma_lognormal(iMode)) )
           
           
        CASE(4)

           nNa = mass_num(:,:,:,isscos)/(1.e-3*wNaCl) ! mol/gridbox
           nCl = nNa
           nSO4= mass_num(:,:,:,iso4cos)/(1.e-3*wSO4)
           nNa2SO4 = MIN(nNa/2.e0, nSO4)
           nNa = nNa - 2.e0*nNa2SO4 ! remaining nNa
           nNaCl = MIN(nCl, nNa)
           nCl = nNaCl ! overshoot nCl is supposed to evaporate
           nH2SO4 = nSO4 - nNa2SO4
           Vol  = ( nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) &
                + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3) &
                + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3) &
                + mass_num(:,:,:,ipomcos)/pom_density &
                + mass_num(:,:,:,isoacos)/soa_density &
                  + mass_num(:,:,:,ibccos)/carbon_density &
                  + mass_num(:,:,:,iducos)/dust_density) ! m3/gridbox
           cmr(iMode,:,:,:) = ( Vol/mass_num(:,:,:,icos_n)*0.75/pi)**(1./3.)*ram2cmr(iMode)! m

           volsu= volsu + nH2SO4*1.e-3*wH2SO4/(dH2SO4*1.e3)
           volna2so4= volna2so4 + nNa2SO4*1.e-3*wNa2SO4/(dNa2SO4*1.e3)
           volpom = volpom + mass_num(:,:,:,ipomcos)/pom_density 
           volsoa = volsoa + mass_num(:,:,:,isoacos)/soa_density 
           volbc = volbc + mass_num(:,:,:,ibccos)/carbon_density 
           volss =volss + nNaCl*1.e-3*wNaCl/(dNaCl*1.e3) 
               
           voldu = voldu + mass_num(:,:,:,iducos)/dust_density

           N_r(i_rad,:,:,:) = N_r(i_rad,:,:,:)+mass_num(:,:,:,mode_tracers(0,iMode)) * 0.5*erfc( log(rc(i_rad)/cmr(iMode,:,:,:)) / &
                sqrt(2.e0)/log(sigma_lognormal(iMode)) )
           
        END SELECT

        volsum=volsum+vol

     end do
     Kappa_r(i_rad,:,:,:) =  (&
           Kap_su*volsu&
          + Kap_na2so4*volna2so4&
          + Kap_pom*volpom&
          + Kap_soa*volsoa&
          + Kap_bc*volbc&
          + Kap_ss*volss&
          + Kap_no3*volno3&
          + Kap_du*voldu&
          ) / volsum
     !where (Kappa_r>1.0)Kappa_r=1.0     
     !  where(Kappa_r < 0.04) Kappa_r = 0.04
    ! Kappa_r(i,:,:,:) =kappa(i,:,:,:)+ Kappa/volsum
        
     
  end do
  deallocate(Vol)
  deallocate(Volsu)
  deallocate(Volss)
  deallocate(Voldu)
  deallocate(Volbc)
  deallocate(Volsoa)
  deallocate(Volpom)
  deallocate(Volna2so4)
  deallocate(Volmsa)
  deallocate(Volno3)
  deallocate(Volsum)

  deallocate(over_rad_frac)

  deallocate(Kappa)
  deallocate(A_Koehler)
  deallocate(N_act)
  deallocate(rc)
  deallocate(nNa)
  deallocate(nCl)
  deallocate(nSO4)
  deallocate(nNaCl)
  deallocate(nNa2SO4)
  deallocate(nH2SO4)

end subroutine Kappa_diag
end module user_output_general