ecearth3/sources/tm5mp/proj/cb05/user_output_aerchemmip.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_AERCHEMMIP
!
! !DESCRIPTION: 
!   isoanus,isoaais,isoaacs,isoacos,isoaaii
!    This module provides the code needed to produce the CMIP6 AERCHEMMIP
!    output from TM5. Code is based on the user_output_aerocom.
!
!    output_aerchemmip_init:
!     - initialise the list of variables (allvars)
!     - initialise the data holder within each variable (2Ddata/3Ddata,...)
!     - initialise the output netcdf files, one for eacht variable
!    output_aerchemmip_accumulate:
!     - do accumulation for all variables and save data to either
!       2Ddata or 3Ddata holder of the variable (excluding optical vars)
!    output_aerchemmip_write
!     - write the monthly variable data to netcdf-file
!     - initialise the dataholders for accumulation of new month
!    output_aerchemmip_write_hourly
!     - write the hourly variable data to netcdf-file
!     - initialise the dataholders for accumulation of new hour
!    output_aerchemmip_write_daily
!     - write the daily variable data to netcdf-file
!     - initialise the dataholders for accumulation of new day
!    write_var
!    
!    calculate_optics
!     - calculate the abss,aods and extinctions
!     - directly accumulate the data containers of these variables
!
!    mode_fraction
!     - calculate the fraction of a M7 mode that is under a size limit

!\\
!\\
! !INTERFACE: 
!
MODULE USER_OUTPUT_AERCHEMMIP
  !
  ! !USES:
  !
  use go,           only : gol, goErr, goPr, goLabel
  use GO,             ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start
  use dims,         only : nregions !=1, support for zooming with larger values not available for AERCHEMMIP
  use optics,       only : wavelendep
  use meteodata,    only : global_lli, levi
  !use chem_param,   only : nmod,xmc2h6,xmc3h6,xmc3h8,xmch4,xmco,xmdms,xmh2o,xmhno3,xmisop,xmno,xmno2,xmo3,xmoh,xmpan,xmso2
  use MDF
  use TM5_DISTGRID, only : dgrid, Get_DistGrid, update_halo, update_halo_iband
  use chem_param
  use m7_data,     only : h2o_mode
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: output_aerchemmip_init
  !public :: output_aerchemmip_step
  public :: output_aerchemmip_write
  public :: output_aerchemmip_write_hourly
  public :: output_aerchemmip_write_6hourly
  public :: output_aerchemmip_write_daily
  public :: output_aerchemmip_done
  public :: accumulate_data
  public :: wdep_out
  logical, public :: aerchemmip_1h = .true.  ! signal for hourly AERCHEMMIP diagnostics
  
  character(len=*), parameter ::  mname = 'user_output_aerchemmip'
  ! max indices, at maximum 7, one for each mode
  integer,parameter :: n_indices=11

  type varfile
     integer                         :: itm5             !
     character(len=16)               :: vname            !
     character(len=64)               :: lname            !
     character(len=11)               :: unit             !
     character(len=10)               :: positive         !
     character(len=130)              :: standard_name    !
     !     real,dimension (:,:),pointer    :: dataZonal           !
     real,dimension (:,:),pointer    :: data2D           !
     real,dimension (:,:,:),pointer  :: data3D           !
     real,dimension (:,:,:),pointer  :: budgetdata          !

     integer                         :: varid            !
     integer                         :: fileunit         ! file unit number 
     character(len=200)              :: filename         !
     integer                         :: dimensions       ! 
     integer                         :: lon_varid       ! 
     integer                         :: lat_varid       ! 
     integer                         :: lev_varid       ! 
     integer                         :: time_varid
     integer                         :: hyam_varid
     integer                         :: hybm_varid
     integer                         :: hyai_varid
     integer                         :: hybi_varid
     integer                         :: bounds_varid
     integer                         :: dims
     character(len=10)               :: freq
     character(len=9)                :: class            ! which class of variable :emi, ddep, wdep,conc,aod,met,crescendo
     integer,dimension(n_indices)         :: indices
     real                            :: xmgas
     character(len=20)               :: table_id
  end type varfile

  type dimdata
     integer               :: nlon             ! size of x dimension of requested field
     integer               :: nlat             ! size of y dimension of requested field
     integer               :: nlev             ! size of 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 dimdata
  type(dimdata)::dimension_data

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

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

  !  type experimentdata
  !     character(len=16)               :: 

  !  end type experimentdata

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

!!!!
  integer::test_fileunit
!!!!

  integer :: n_vars=0
  type(varfile), dimension(:), allocatable :: allvars
  type(varfile), dimension(:), allocatable :: fixedvars
  integer :: n_var_max=300
  integer :: n_fixed=3
  integer, public :: n_days_in_month
  real :: sfo3_molmol
  character(len=20), public :: aerchemmip_exper ! AeroCom experiment name

  integer, save             :: od550aer, &
       ec550aer,&
       ec550aermon,&
       od550aerdaily, &
       od550so4, &
       od550bc,   &
       od550oa,&
       od550soa,&
       od550ss,&
       od550dust,&
       od440dustday,&
       od550dust2dday,&
       od550dust3dday,&
       od550no3,&
       od550aerh2o,&
       od550lt1aer,&
       abs550aer,&
       od440aer,&
       od350aer,&
       od870aer,&
       od440aerhr,&
       od440aermonthly,&
       od440aerdaily,&
       od550aerhr,&
       areacella,&
       sftlf,&
       orog

  integer                                   ::  fid ! file id for IF_NOTOK_MDF macro
  integer                                   :: access_mode
  integer                                   :: accumulation_mon,accumulation_day,accumulation_hr,accumulation_6hr
  integer                                   :: timeidx_mon,timeidx_hr,timeidx_day,timeidx_6hr
  logical,public::crescendo_out=.false.

  integer,parameter::iemiunit=1
  integer,parameter::iddepunit=1 !same dimensions as emi
  integer,parameter::iwdepunit=1 !same dimensions as emi
  integer,parameter::iprod3dunit=2
  integer,parameter::immrunit=3
  integer,parameter::idimensionlessunit=4
  integer,parameter::iheightunit=5
  integer,parameter::itempunit=6
  integer,parameter::io3unit=7
  integer,parameter::ipresunit=8
  integer,parameter::ivmrunit=9
  integer,parameter::irateunit=10
  integer,parameter::iloadunit=11
  integer,parameter::iextunit=12
  integer,parameter::iccunit=13
  integer,parameter::im3unit=14
  integer,parameter::imassunit=15

  character(len=11),dimension(15),parameter::units=(/'kg m-2 s-1','kg m-3 s-1','kg kg-1','1','m','K','DU','Pa','mole mole-1',&
       's-1','kg m-2','m-1','cm-3','m-3','kg'/)  

  !CRESCENDO
  !3D

  Character(len=11),dimension(40),parameter   :: crescendo3D=    (/'nd1', 'nd2', 'nd3', 'nd4', 'nd5', 'nd6', 'nd7', &
       'mmrsu1', 'mmrsu2', 'mmrsu3', 'mmrsu4', 'mmrsoa1', 'mmrsoa2', 'mmrsoa3', 'mmrsoa4', 'mmrsoa5', 'mmroa2', &
       'mmroa3', 'mmroa4', 'mmroa5', 'mmrbc2', 'mmrbc3', 'mmrbc4', 'mmrbc5', 'mmrss3', 'mmrss4', 'mmrdu3', 'mmrdu4', &
       'mmrdu6', 'mmrdu7', 'mmraerh2o_1', 'mmraerh2o_2', 'mmraerh2o_3', 'mmraerh2o_4', 'mono', 'nh3', 'ndtot', &
       'ec550aer', 'chepsoa3d','emilnox'/)

  Character(len=11),dimension(40),parameter   :: crescendo3Dunit=(/units(im3unit), units(im3unit), units(im3unit), &
       units(im3unit), units(im3unit), units(im3unit), units(im3unit), units(immrunit), units(immrunit), units(immrunit), &
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),&
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), &
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), &
       units(immrunit), units(immrunit), units(immrunit), units(ivmrunit), units(ivmrunit), units(im3unit), units(iextunit), &
       units(iemiunit),units(iemiunit)/)

  Character(len=11),dimension(2),parameter   :: crescendo3Dday=(/'co', 'od5503ddust'/)
  Character(len=11),dimension(2),parameter   :: crescendo3Ddayunit=(/units(ivmrunit),units(idimensionlessunit)/)


  !2D
  !mon
  Character(len=8),dimension(14),parameter   :: crescendo2Dmon=(/'drydms', 'wetdms', 'wetno3', 'dryhno3', 'wethno3', &
       'dryno2', 'dryno', 'drypan', 'emimono', 'dmsos', 'seddust', 'uas', 'vas', 'sfcWind'/)
  Character(len=11),dimension(14),parameter   :: crescendo2Dmonunit=(/units(iddepunit), units(iddepunit), units(iddepunit), &
       units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), &
       'kg m-3',units(iddepunit),'m s-1', 'm s-1', 'm s-1'/)


  ! 2d
  ! 6hr
  Character(len=16),dimension(19),parameter   :: crescendo2D6hr=(/'sfmmrso4', 'sfmmrss', 'sfmmroa', 'sfmmrsoa', 'sfmmrbc', &
       'sfmmrdustabv1', 'sfmmrdustabv10', 'sfmmrdustbel1', 'sfdms', 'sfisop', 'sfmono', 'emidms', 'emiss', 'emioa', &
       'emiisop', 'emimono', 'sfndtot', 'sfnd100', 'chepsoa2d'/)
  Character(len=11),dimension(19),parameter   :: crescendo2D6hrunit=(/units(immrunit), units(immrunit), units(immrunit),&
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(ivmrunit), &
       units(ivmrunit), units(ivmrunit), units(iemiunit), units(iemiunit), units(iemiunit), units(iemiunit), &
       units(iemiunit), units(im3unit), units(im3unit),'kgm-2s-1'/)
  !2d
  !day
  Character(len=16),dimension(52),parameter   :: crescendo2Dday=(/'sfnd1', 'sfnd2', 'sfnd3', 'sfnd4', 'sfnd5', 'sfnd6', &
       'sfnd7', 'sfmmrsu1', 'sfmmrsu2', 'sfmmrsu3', 'sfmmrsu4', 'sfmmrsoa1', 'sfmmrsoa2', 'sfmmrsoa3', 'sfmmrsoa4', &
       'sfmmrsoa5', 'sfmmroa2', 'sfmmroa3', 'sfmmroa4', 'sfmmroa5', 'sfmmrbc2', 'sfmmrbc3', 'sfmmrbc4', 'sfmmrbc5', &
       'sfmmrss3', 'sfmmrss4', 'sfmmrdu3', 'sfmmrdu4', 'sfmmrdu6', 'sfmmrdu7', 'sfmmraerh2o_1', 'sfmmraerh2o_2', &
       'sfmmraerh2o_3', 'sfmmraerh2o_4', 'sfmmrno3', 'sfmmrnh4', 'sfmmrnh3', 'sfsh', 'od440aer', 'od440dust', &
       'od550dust', 'depdustbel1', 'depdustabv1', 'depdustabv10', 'sfmmrdust', 'drynh4', 'wetnh4', 'dryno3', &
       'wetno3', 'dryhno3','drydust','wetdust'/) 

  Character(len=11),dimension(52),parameter   :: crescendo2Ddayunit=(/units(im3unit), units(im3unit), units(im3unit), &
       units(im3unit), units(im3unit), units(im3unit), units(im3unit), units(immrunit), units(immrunit), units(immrunit), &
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),&
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), &
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),&
       units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit), units(immrunit),&
       units(idimensionlessunit), units(idimensionlessunit), units(idimensionlessunit), units(idimensionlessunit), units(iddepunit),& 
       units(iddepunit), units(iddepunit), units(immrunit), units(iddepunit), units(iddepunit), units(iddepunit), units(iddepunit), &
       units(iddepunit), units(iddepunit), units(iddepunit)/)

  !2d
  !hr
  Character(len=9), dimension(5),  parameter   :: crescendo2Dhr=(/'od550aer', 'od440aer', 'sfno', 'sfnh3', 'sfhno3'/)

  Character(len=11), dimension(5), parameter   :: crescendo2Dhrunit=(/units(idimensionlessunit), units(idimensionlessunit), &
       units(ivmrunit), units(ivmrunit), units(ivmrunit)/)

  !Character(len=11),dimension(6,2),parameter :: crescendo2Dtest=reshape(&
  !    (/'od550aer', 'od440aer', 'sfno'         ,'sfnh3'        ,'sfhno3'       ,'chepsoa2d' ,&
  !      '1'       ,'1'       ,units(ivmrunit), units(ivmrunit), units(ivmrunit),'kgm-2s-1'/),(/6,2/))


  !METEO
  !3D
  Character(len=7),dimension(11),parameter   :: meteo3D=(/'ta', 'pfull', 'phalf', 'hus', 'zg', 'airmass', 'emilnox', &
       'jno2', 'ua', 'va', 'wa'/)
  Character(len=11),dimension(11),parameter   :: meteo3Dunit=(/units(itempunit), units(ipresunit), units(ipresunit), &
       units(idimensionlessunit), units(iheightunit), units(iloadunit), 'kg N/m2/month', units(irateunit),'ms-1', 'ms-1', 'ms-1'/)
  !2D
  Character(len=7),dimension(9),parameter   :: meteo2D=(/'ps', 'ptp', 'tatp', 'ztp', 'bldep', 'pr', 'tropoz', 'toz', 'albsrfc'/)
  Character(len=11),dimension(9),parameter   :: meteo2Dunit=(/units(ipresunit), units(ipresunit), units(itempunit), &
       units(iheightunit), units(iheightunit), units(iemiunit), units(io3unit), units(io3unit), units(idimensionlessunit)/)

  !OPTICS
  Character(len=11),dimension(13),parameter   :: opticscomp=(/'od550aer', 'od550so4', 'od550bc', 'od550oa', 'od550soa',&
       'od550ss', 'od550dust', 'od550no3', 'od550aerh2o', 'od550lt1aer', 'od440aer', 'od870aer', 'abs550aer'/)!
  !AEROSOL compounds
  Character(len=3),dimension(6),parameter   :: comp=(/'BC', 'POM', 'SO4', 'DU', 'SS', 'SOA'/)!SOA

  !CHEMICAL
  Character(len=6),dimension(2),parameter  :: checomp=(/'aqpso4', 'gpso4'/)
  Character(len=6),dimension(1),parameter  :: chepcomp=(/'soa'/)
  Character(len=7),dimension(4),parameter  :: o3chepcomp=(/'o3loss', 'o3prod','lossch4','lossco'/)

  !Emon
  Character(len=9),dimension(8),parameter  :: emonout=(/'flashrate', 'depdust','seddustCI','md','loaddust','concdust','conccn','sconcss'/)
  Character(len=11),dimension(8),parameter  :: emonoutunit=(/'km-2 s-1', units(iddepunit),units(iddepunit),'m','kg m-2','kg m-3','m-3','kg m-3'/)

  !BUDGET (EMI+REMOVAL)
  Character(len=4),dimension(14),parameter  :: emicomp=(/'bvoc', 'co', 'dms', 'isop', 'nox', 'nh3', 'oa', 'so2', 'bc', &
       'so4', 'dust', 'ss', 'voc','poa'/)
  Character(len=4),dimension(12),parameter  :: ddepcomp=(/'nh3', 'noy', 'o3', 'oa', 'so2', 'bc', 'so4', 'dust', 'ss', &
       'soa', 'nh4', 'no3'/)
  Character(len=4),dimension(10),parameter   :: wdepcomp=(/'bc', 'dust', 'nh3', 'nh4', 'noy', 'oa', 'so2', 'so4', 'ss', 'soa'/)

  !MIXING ratios
  Character(len=6),dimension(13),parameter   :: aerommrcomp=(/'bc', 'dust', 'nh3', 'nh4', 'no3', 'oa', 'so4', 'ss', 'pm1', &
       'pm2p5', 'pm10', 'aerh2o', 'soa'/)
  Character(len=8),dimension(20),parameter   :: gascomp=(/'c2h6', 'c3h6', 'c3h8', 'ch3coch3', 'ch4', 'co', 'co2', 'dms', &
       'h2o', 'hno3', 'isop', 'no', 'no2', 'o3', 'oh', 'pan', 'so2', 'voc', 'hcho', 'o3ste'/)  
  !Molecular weights
  real,dimension(20),parameter   :: xmgascomp=(/xmc2h6, xmc3h6, xmc3h8, xmacet, xmch4, xmco, -1.0, xmdms, xmh2o, xmhno3, &
       xmisop, xmno, xmno2, xmo3, xmoh, xmpan, xmso2, 1.0, xmch2o,xmo3/)!VOC=1, 
  !hcho=ch2o in TM5, but output for hcho is needed.

  ! HOURLY and 6-HOURLY
  character(len=8),  dimension(2), parameter :: hr6_var=(/'ps','ec550aer'/)
  character(len=11), dimension(2), parameter :: hr6_unit=(/units(ipresunit), units(iextunit)/)
  character(len=8),  dimension(5), parameter :: hourly_var=(/'ps', 'sfno2', 'sfo3', 'sfpm25', 'tas'/)
  character(len=11), dimension(5), parameter :: hourly_varunit=(/units(ipresunit), units(ivmrunit), units(ivmrunit), &
       units(immrunit), units(itempunit)/)

  !DAILY
  character(len=8),dimension(10),parameter:: daily_var=(/'od550aer', 'toz', 'maxpblz', 'minpblz', 'tasmin', 'tasmax', &
       'pr', 'sfo3max', 'tas','ps'/)
  character(len=11),dimension(10),parameter:: daily_varunit=(/units(idimensionlessunit), units(io3unit), &
       units(iheightunit), units(iheightunit), units(itempunit), units(itempunit), units(iemiunit), units(ivmrunit), units(itempunit),units(ipresunit)/)

  !ZONAL
  character(len=6),dimension(8),parameter:: zonal_var=(/'ch4', 'hno3', 'ho2', 'noy', 'ta', 'zg', 'o3', 'oh'/)
  character(len=11),dimension(8),parameter:: zonal_varunit=(/units(ivmrunit), units(ivmrunit), units(ivmrunit), &
       units(ivmrunit), units(itempunit), units(iheightunit), units(ivmrunit), units(ivmrunit)/)
  integer,dimension(8),parameter:: zonal_idx=(/ich4,ihno3,iho2,-1,-1,-1,io3,ioh/)


  !AERCHEMMIP global attributes that might change with run or something else
  character(len=3),parameter::grid='3x2'!'250 km'
  character(len=3),parameter::grid_label='gn'!'gnz' for zonal means
  character(len=300),parameter::source='EC-Earth3-AerChem (2017): atmosphere: IFS cy36r4 (TL255, linearly &
       &reduced Gaussian grid equivalent to 512 x 256, 91 levels, top level: 0.01 hPa);atmospheric_chemistry: &
       &TM5 (3 deg. (long.) x 2 deg. (lat.), 34 levels, top level: 0.1 hPa; aerosol: TM5'
  character(len=17),parameter::source_id='EC-Earth3-AerChem'
  character(len=20),public ::source_type!='AOGCM CHEM AER' !or 'AGCM CHEM AER' for amip-runs
  character(len=60),public ::realm
  character(len=60),public::experiment_id
  character(len=60),public::experiment
  character(len=1),public::realization_i='1'
  character(len=1),public::physics_i='1'
  character(len=1),public::forcing_i='1'
  character(len=1),public::initialization_i='1'
  integer, public:: aerchemmip_dhour
  ! Timers
  integer :: itim_init, itim_addvar, itim_write, itim_accu, itim_attr, itim_accu_opt, itim_write_hour, itim_write_day, &
       itim_write_mon, itim_write_gather
contains

  subroutine output_aerchemmip_init(status)
    ! Open files
    ! allocate dataholders
    use dims,            only : newsrun,itau,mlen
    use global_data,     only : outdir
    use datetime,        only : tau2date, date2tau
    use partools,        only : MPI_INFO_NULL, localComm
    use optics,          only : Optics_Init
    use sedimentation,   only : ised,nsed
    use partools    ,    only : isRoot,myid
    use global_data,     only : region_dat
    use tm5_distgrid,    only : gather
    use meteodata      , only : lsmask_dat,oro_dat
    use Binas      ,     only : grav

    implicit none

    ! OUTPUT parameters
    integer, intent(out)     :: status

    ! LOCAL parameters
    integer  :: region !iterator for regions
    integer  :: nlon_region
    integer  :: nlat_region
    integer  :: nlev_region ! also global

    integer  :: j_var
    integer  :: itrac
    integer  :: i_sed
    integer  :: i,i1,i2,j1,j2,k,j,imr,jmr
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_init'
    ! FIXED  VARS
    real, dimension(:),pointer :: dxyp 
    real, allocatable :: arr2d(:,:)
    real ::xmcomp

    call goLabel(rname)

    ! define timers:
    call GO_Timer_Def( itim_init, 'output aerchemmip init', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write, 'output aerchemmip write', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_gather, 'output aerchemmip write gather', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_day, 'output aerchemmip write day', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_hour, 'output aerchemmip write hour', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_write_mon, 'output aerchemmip write mon', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_accu, 'output aerchemmip accu', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_accu_opt, 'output aerchemmip accu _optics', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_attr, 'output aerchemmip attr', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_addvar, 'output aerchemmip addvar', status )
    IF_NOTOK_RETURN(status=1)

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

    if (newsrun) then
      !optics?
      ! 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)

      !TB Have to keep insitu part, since optics-module uses it for comparisons.

      allocate( wdep_out( 3 ) )
      wdep_out(1)%wl = 0.550 ; wdep_out(1)%split = .true.  ; wdep_out(1)%insitu = .false. 
      wdep_out(2)%wl = 0.440 ; wdep_out(2)%split = .true. ; wdep_out(2)%insitu = .false. 
      wdep_out(3)%wl = 0.870 ; wdep_out(3)%split = .false. ; wdep_out(3)%insitu = .false.
      !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)
      accumulation_mon=0.0
      accumulation_hr=0.0
      accumulation_6hr=0.0
      accumulation_day=0.0
      region=1
      
      ! intermediate structures for budgets
      allocate(drydepos(region)%f2dslast(i1:i2,j1:j2,18))
      allocate(wetdepos(region)%f2dslast(i1:i2,j1:j2,13))
      allocate(emission(region)%f2dslast(i1:i2,j1:j2,13))
      
      ! these here are the initial budgets (monthly): 0.0
      drydepos(region)%f2dslast = 0.0
      wetdepos(region)%f2dslast = 0.0
      emission(region)%f2dslast = 0.0
      imr = global_lli(1)%nlon
      jmr = global_lli(1)%nlat

      ! for areacella, orog, sftlf
      if (isRoot) then
        allocate( arr2d(imr,jmr)       )
      else
        allocate( arr2d(1,1)   )
      endif
      arr2d(:,:)=0.0

      ! for monthly output, initialise with 31 for january
      n_days_in_month=31

    end if

    call GO_Timer_Start( itim_init, status )
    IF_NOTOK_RETURN(status=1)
    ! AERCHEMMIP only available for global-> region=1
    region=1

    ! Initialize grid definitions
    nlon_region = global_lli(region)%nlon
    nlat_region = global_lli(region)%nlat
    nlev_region = levi%nlev
    dimension_data%nlon= nlon_region
    dimension_data%nlat= nlat_region
    dimension_data%nlev= nlev_region
    allocate(dimension_data%lon(nlon_region))       
    allocate(dimension_data%lat(nlat_region))       
    allocate(dimension_data%lev(nlev_region))       
    dimension_data%lon=global_lli(region)%lon_deg
    dimension_data%lat=global_lli(region)%lat_deg

    ! initialise output timeidx used for keeping track which output steps is written
    timeidx_mon=1
    timeidx_day=1
    timeidx_hr=1
    timeidx_6hr=1
    ! allocate room for variables
    allocate(allvars(n_var_max))       
    allocate(fixedvars(n_fixed))       

    if (crescendo_out)then

      do i=1,size(crescendo3D)
        if (trim(crescendo3D(i))=='mono')then
          xmcomp=xmterp
        else if  (trim(crescendo3D(i))=='nh3')then
          xmcomp=xmnh3
        else
          write(gol,*) 'CRESCENDO 3D monthly with negative molar mass'             
          xmcomp=-1.0
        end if
        call add_variable(-1,trim(crescendo3D(i)),trim(crescendo3D(i)),crescendo3Dunit(i),3,status,'crescendo','AERmon',xmcomp)
      end do

      do i=1,size(crescendo3Dday)
        if (trim(crescendo3Dday(i))=='co')then
          xmcomp=xmco
        else
          write(gol,*) 'CRESCENDO 3D daily with negative molar mass'             
          xmcomp=-1.0
        end if
        call add_variable(-1,trim(crescendo3Dday(i)),trim(crescendo3Dday(i)),crescendo3Ddayunit(i),3,status,'crescendo','AERday',xmcomp)
      end do
      
      do i=1,size(crescendo2D6hr)
        if (trim(crescendo2D6hr(i))=='sfdms')then
          xmcomp=xmdms
        else if (trim(crescendo2D6hr(i))=='sfisop')then
          xmcomp=xmisop
        else if (trim(crescendo2D6hr(i))=='sfmono')then
          xmcomp=xmterp
        else
          write(gol,*) 'CRESCENDO 2D 6hr with negative molar mass'             
          write(gol,*) 'gascomp with negative molar mass'
          xmcomp=-1.0
        end if
        call add_variable(-1,trim(crescendo2D6hr(i)),trim(crescendo2D6hr(i)),crescendo2D6hrunit(i),2,status,'crescendo','AER6hr',xmcomp)
      end do

      do i=1,size(crescendo2Dmon)
        call add_variable(-1,trim(crescendo2Dmon(i)),trim(crescendo2Dmon(i)),crescendo2Dmonunit(i),2,status,'crescendo','AERmon',-1.0)
      end do

      do i=1,size(crescendo2Dhr)
        if (trim(crescendo2Dhr(i))=='sfno')then
          xmcomp=xmno
        else if (trim(crescendo2Dhr(i))=='sfnh3')then
          xmcomp=xmnh3
        else if (trim(crescendo2Dhr(i))=='sfhno3')then
          xmcomp=xmhno3
        else
          ! -1 so that missing molar mass will be noticed easily
          write(gol,*) 'CRESCENDO 2D hr with negative molar mass'             
          xmcomp=-1.0
        end if
        call add_variable(-1,trim(crescendo2Dhr(i)),trim(crescendo2Dhr(i)),crescendo2Dhrunit(i),2,status,'crescendo','AERhr',xmcomp)
      end do

      do i=1,size(crescendo2Dday)
        call add_variable(-1,trim(crescendo2Dday(i)),trim(crescendo2Dday(i)),crescendo2Ddayunit(i),2,status,'crescendo','AERday',-1.0)
        !call add_variable(-1,trim(crescendo2Dday_new(i,1)),trim(crescendo2Dday_new(i,1)),crescendo2Dday_new(i,2),2,status,'crescendo','AERday',-1.0)
      end do
    end if

    do i=1,size(hr6_var)
      if (trim(hr6_var(i))=='ec550aer')then
        call add_variable(-1,trim(hr6_var(i)),'optics '//trim(hr6_var(i)), hr6_unit(i),3,status,'optics','AER6hr',-1.0)
      else
        call add_variable(-1,trim(hr6_var(i)),trim(hr6_var(i)),hr6_unit(i),2,status,'ps6h','AER6hr',-1.0)
      endif
    end do

    ! add deposition variables
    do i=1,size(emicomp)
      call add_variable(-1,'emi'//trim(emicomp(i)),'emission '//trim(emicomp(i)), units(iemiunit),2,status,'emi','AERmon',-1.0)
    end do

    ! add 3D chemical production fields 
    do i=1,size(checomp)
      call add_variable(-1,'che'//trim(checomp(i)),'chemical production of '//trim(checomp(i)), units(iemiunit),3,status,'chep','AERmon',-1.0)
    end do

    ! add 3D ozone prod loss
    do i=1,size(o3chepcomp)
      call add_variable(-1,trim(o3chepcomp(i)),'tendency_'//trim(o3chepcomp(i)),'mol m-3 s-1',3,status,'chep','AERmon',-1.0)
    end do

    ! add 2D chemical production fields 
    do i=1,size(chepcomp)
      call add_variable(-1,'chep'//trim(chepcomp(i)),'chemical production of '//trim(chepcomp(i)), units(iemiunit),2,status,'chep','AERmon',-1.0)
    end do

    ! add Emon fields
    do i=1,size(emonout)
      select case (trim(emonout(i)))
      case( 'md','concdust','conccn')
        call add_variable(-1,trim(emonout(i)),trim(emonout(i)), emonoutunit(i),3,status,'emon','Emon',-1.0)
      case('flashrate','depdust','seddustCI','loaddust','sconcss')
        call add_variable(-1,trim(emonout(i)),trim(emonout(i)), emonoutunit(i),2,status,'emon','Emon',-1.0)
      end select
    end do

    ! add dry deposition fields
    do i=1,size(ddepcomp)
      call add_variable(-1,'dry'//trim(ddepcomp(i)),'dry deposition '//trim(ddepcomp(i)), units(iddepunit),2,status,'dry','AERmon',-1.0)
    end do

    ! add wetdep variables
    do i=1,size(wdepcomp)
      call add_variable(-1,'wet'//trim(wdepcomp(i)),'wet deposition '//trim(wdepcomp(i)), units(iwdepunit),2,status,'wet','AERmon',-1.0)
    end do

    ! add optics fields
    do i=1,size(opticscomp)
      !if (trim(opticscomp(i))=='ec550aer') then
      !   call add_variable(-1,trim(opticscomp(i)),'optics '//trim(opticscomp(i)), units(iextunit),3,status,'optics','AER6hr',-1.0)
      !else
      call add_variable(-1,trim(opticscomp(i)),'optics '//trim(opticscomp(i)), units(idimensionlessunit),2,status,'optics','AERmon',-1.0)
    end do

    ! Aerosol species mass mixing ratios
    do i=1,size(aerommrcomp)
      call add_variable(-1,'mmr'//trim(aerommrcomp(i)),'mass mixing ratio of '//trim(aerommrcomp(i)), units(immrunit),3,status,'mmr','AERmon',-1.0)
    end do

    ! Gas-phase species volume mixingratios
    do i=1,size(gascomp)
      if (xmgascomp(i)>0.0) then
        call add_variable(-1,trim(gascomp(i)),'volume mixing ratio of '//trim(gascomp(i)), units(ivmrunit),3,status,'gas','AERmon',xmgascomp(i))
      else
        write(gol,*) 'gascomp with negative molar mass'
      end if
    end do

    ! add meterorological fields
    do i=1,size(meteo3D)
      call add_variable(-1,trim(meteo3D(i)),trim(meteo3D(i)),meteo3Dunit(i),3,status,'meteo3D','AERmon',-1.0)
    end do

    ! surface meteorological fields
    do i=1,size(meteo2D)
      call add_variable(-1,trim(meteo2D(i)),trim(meteo2D(i)),meteo2Dunit(i),2,status,'meteo2D','AERmon',-1.0)
    end do

    ! add hourly output
    if ( aerchemmip_1h ) then
      do i=1,size(hourly_var)
        call add_variable(-1,trim(hourly_var(i)),trim(hourly_var(i)),hourly_varunit(i),2,status,'sf1h','AERhr',-1.0)
      end do      
    end if

    ! add daily fields
    do i=1,size(daily_var)
      call add_variable(-1,trim(daily_var(i)),trim(daily_var(i)),daily_varunit(i),2,status,'sf1d','AERday',-1.0)
    end do

    ! add zonal fields, some fields are repeated (3d/zonal)
    do i=1,size(zonal_var)
      call add_variable(zonal_idx(i),trim(zonal_var(i)),trim(zonal_var(i)),zonal_varunit(i),3,status,'zonal','AERmonZ',-1.0)
    end do

    ! Fixed fields
    call add_variable(-1,'areacella','cell area','m2',2,status,'fixed','AERfx',-1.0)
    call add_variable(-1,'orog','surface_altitude','m',2,status,'fixed','AERfx',-1.0)
    call add_variable(-1,'sftlf','land_area_fraction','1',2,status,'fixed','AERfx',-1.0)

    ! Initialize a single file for each variables as per AERCHEMMIP specification
    do j_var = 1, n_vars

      ! overwrite existing files (clobber)
      if (isroot)call MDF_Create( allvars(j_var)%filename, MDF_NETCDF4, MDF_REPLACE, allvars(j_var)%fileunit, status )
      IF_NOTOK_RETURN(status=1)

      ! write grid dimension attributes
      if (isroot) call write_dimensions(status, j_var)
      IF_NOTOK_RETURN(status=1)

      ! write global attributes
      if (isroot) call write_attributes(status, j_var)
      IF_NOTOK_RETURN(status=1)

      ! write dimension variables
      if (isroot) call write_var(status,j_var)
      IF_NOTOK_RETURN(status=1)

      ! Fixed fields
      if (allvars(j_var)%table_id=='AERfx')then

        if (trim(allvars(j_var)%vname)=='areacella')then
          ! Gridbox area
          dxyp   => region_dat(1)%dxyp
          do j=j1,j2
            allvars(j_var)%data2D(i1:i2,j)=dxyp(j)
          end do
          call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        else  if (trim(allvars(j_var)%vname)=='orog')then
          ! Orography
          allvars(j_var)%data2D(i1:i2,j1:j2) =oro_dat(region)%data(i1:i2,j1:j2,1)/grav
          call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        else  if (trim(allvars(j_var)%vname)=='sftlf')then
          ! LSM
          allvars(j_var)%data2D(i1:i2,j1:j2)=lsmask_dat(1)%data(i1:i2,j1:j2,1)/100.
          call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1/), count=(/imr,jmr/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        end if
      end if
    end do

    deallocate(arr2d)
    call GO_Timer_End( itim_init, status )
    IF_NOTOK_RETURN(status=1)       
    call goLabel()
    status = 0

  end subroutine output_aerchemmip_init

  subroutine output_aerchemmip_write(region,newhour,status)
    use MeteoData,      only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,&
         blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat   !
    use global_data,    only : conv_dat
    use GO,             only : TDate, NewDate
    use go_date,only: days_in_month!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use ebischeme,      only : AC_diag_prod,AC_O3_lp,AC_loss
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    use partools    , only : isRoot,myid

    !    use domain_decomp, only: gather
    implicit none
    logical,intent(in) ::newhour
    integer,intent(out)::status
    integer::region
    integer:: j_var
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_write'
    real, allocatable ::  arr3d(:,:,:),arr3dh(:,:,:),arr2d(:,:)
    integer, dimension(6) :: curdate    
    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    real                          :: rangemon
    type(Tdate)::curdate_tdate
    call goLabel(rname)

    call GO_Timer_Start( itim_write_mon, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
      write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!'
      call goErr
      status=1
      return
    end if

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

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

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 
    ! returns the difference to current time, beginning of new step
    ! 
    call date2tau( time_reftime6, itauref )

    ! calculate reftime as fractional days from itauref, hence division by 86400
    ! call date2tau( idater, itaucur )
    ! itau is the 1st day of new month at 00:00 so we need to fix the reftime back half a month (15th day)
    ! ((cursecond - reftimesecond) / seconds_in_day) - days_in_last_month + 15days

    !reftime = (itau - itauref -n_days_in_month*24*3600 + 15*24*3600) / 86400.  
    reftime = (itau - itauref ) / 86400.  
    !get current date in integers
    call tau2date(itau, curdate)
    ! create a TDATE date variable of the previous month (curdate(3)-1)

    curdate_tdate=newdate(curdate(1),curdate(2),curdate(3)-1,curdate(4),curdate(5),curdate(6),calender='greg')
    ! get days in month and save for next step
    n_days_in_month=days_in_month(curdate_tdate)
    ! change reftime to beginning of last month (the month data is from)
    reftime=reftime-n_days_in_month
    !length of the month-1s(in days) for the time_bounds
    rangemon=n_days_in_month !-(1.0/86400.0)

    ! allocate 3D and 4D global arrays for gathering data
    ! only root needs the full array, but it must be allocated in all tasks
    if (isRoot) then
      allocate( arr3d(imr,jmr,lmr)       )
      allocate( arr3dh(imr,jmr,lmr+1)       )
      allocate( arr2d(imr,jmr)       )
    else
      allocate( arr3d(1,1,1)   )
      allocate( arr3dh(1,1,1)   )
      allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0
    arr3dh(:,:,:)=0.0
    do j_var =1,n_vars
      ! hourly and daily variables are handled separately
      if (allvars(j_var)%table_id=='AERhr'.or.allvars(j_var)%table_id=='AER6hr'.or.&
           allvars(j_var)%table_id=='AERday'.or.allvars(j_var)%table_id=='AERfx')then

        cycle
      end if
      if (allvars(j_var)%dims==2)then !2D


        if (trim(allvars(j_var)%vname)/='minpblz'.and.trim(allvars(j_var)%vname)/='tasmin'.and. &
             trim(allvars(j_var)%vname)/='maxpblz'.and.trim(allvars(j_var)%vname)/='tasmax')then

          allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_mon)
        end if
        call GO_Timer_Start( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)

        call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0,status)
        call GO_Timer_End( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)

        IF_NOTOK_RETURN(status=1)

        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1,timeidx_mon/), &
             count=(/imr,jmr,1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        ! Zero out the accumulated and written data for the next interval
        if (trim(allvars(j_var)%vname)=='minpblz' .or. trim(allvars(j_var)%vname)=='tasmin') then
          ! put high number so simple comparison is needed for finding minimum
          allvars(j_var)%data2D(i1:i2,j1:j2)=1e10
        else
          allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
        end if
      else !3D
        if (trim( allvars(j_var)%vname)=='phalf') then !lmr+1
          allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr+1)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr+1)/real(accumulation_mon) 
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          call gather( dgrid(1), allvars(j_var)%data3D , arr3dh(:,:,:), 0, status)
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)

          if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3dh , status, start=(/1,1,1,timeidx_mon/), &
               count=(/imr,jmr,lmr+1,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        else
          allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_mon)
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)

          if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3d , status, start=(/1,1,1,timeidx_mon/), &
               count=(/imr,jmr,lmr,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        end if

        ! Zero out the accumulated and written data for the next interval
        allvars(j_var)%data3D(i1:i2,j1:j2,:)=0.0

      end if
      !end if
      ! write the date for timestep
      if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+real(rangemon/2)/) , status, start=(/timeidx_mon/), count=(/1/) )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+rangemon/) , status, &
           start=(/1,timeidx_mon/), count=(/2,1/) )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)



    end do
    deallocate( arr3d,arr3dh,arr2d)
    ! change time index
    timeidx_mon= timeidx_mon + 1
    ! accululated time to zero
    accumulation_mon=0
    ! zero out the chemical production (for mongthly output)
    !AC_diag_prod(region)%prod(i1:i2,j1:j2,:,1:3)=0.0
    ! zero out the chemical production 
    !AC_O3_lp(region)%prod(i1:i2,j1:j2,:,1:2)=0.0
    ! zero out the chemical production 
    !AC_loss(region)%prod(i1:i2,j1:j2,:,1:2)=0.0

    !status=1
    !return
    call GO_Timer_End( itim_write_mon, status )
    IF_NOTOK_RETURN(status=1)


    call goLabel()
    status = 0
  end subroutine output_aerchemmip_write
  
  subroutine output_aerchemmip_write_daily(region,newday,status)
    use MeteoData,      only : temper_dat, sst_dat, albedo_dat, ci_dat, lsp_dat, cp_dat, ssr_dat, str_dat, &
         blh_dat, gph_dat, lwc_dat, iwc_dat, cco_dat, cc_dat, humid_dat, m_dat, phlb_dat, sp_dat   !
    use meteodata   , only : global_lli, levi
    use partools    , only : isRoot,myid
    use GO,             only : TDate, NewDate!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    implicit none
    logical,intent(in) ::newday
    integer,intent(out)::status
    integer::region
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_write'
    integer:: j_var

    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    real                          :: rangeday ! for bounds
    ! root writes netcdf arrays
    real, allocatable :: arr3d(:,:,:), arr2d(:,:)    
    integer:: imr_f,jmr_f,lmr_f

    call goLabel(rname)

    call GO_Timer_Start( itim_write_day, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
      write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!'
      call goErr
      status=1
      return
    end if


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


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

    ! allocate 3D and 4D global arrays for gathering data
    if (isRoot) then
      allocate( arr3d(imr,jmr,lmr)       )
      allocate( arr2d(imr,jmr)       )
    else
      allocate( arr3d(1,1,1)   )
      allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 

    call date2tau( time_reftime6, itauref )
    ! calculate reftime as fractional days from itauref, hence division by 86400
    ! call date2tau( idater, itaucur )
    reftime = (itau - itauref) / 86400. - 1.0
    !23h59 as days
    rangeday=1.0!(23.0*3600.0+59.0*60.0+59.0)/86400.0
    do j_var =1,n_vars
      if (allvars(j_var)%table_id/='AERday')then 
        cycle
      end if
      if (allvars(j_var)%dims==2)then

        if ( trim(allvars(j_var)%vname)/='minpblz' .and. trim(allvars(j_var)%vname)/='tasmin' .and. &
             trim(allvars(j_var)%vname)/='maxpblz' .and. trim(allvars(j_var)%vname)/='tasmax' .and. &
             trim(allvars(j_var)%vname)/='sfo3max' ) then
          allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_day)
        end if
        call GO_Timer_Start( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)
        call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0, status)
        IF_NOTOK_RETURN(status=1)
        call GO_Timer_End( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)

        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr2d, status, start=(/1,1,timeidx_day/), count=(/imr,jmr,1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        if (trim(allvars(j_var)%vname)=='minpblz' .or. trim(allvars(j_var)%vname)=='tasmin') then
          ! put high number so simple comparison is needed for finding minimum
          allvars(j_var)%data2D(i1:i2,j1:j2)=1e10
        else
          ! Zero out the accumulated and written data for the next interval
          allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
        end if
      else

        allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_day)
        !end if
        call GO_Timer_Start( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)
        call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
        call GO_Timer_End( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)          !if (trim(allvars(j_var)%vname)=='od5503ddust')then
        IF_NOTOK_RETURN(status=1)
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid, arr3d, status, start=(/1,1,1,timeidx_day/), count=(/imr,jmr,lmr,1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)


        allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0

      end if


      ! write the date for timestep
      if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+0.5/) , status, start=(/timeidx_day/), count=(/1/) )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+ rangeday/) , status, start=(/1,timeidx_day/), count=(/2,1/) )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

    end do
    deallocate(arr3d, arr2d)
    ! Timeindex to next day
    timeidx_day= timeidx_day + 1
    ! daily accumulated time to zero
    accumulation_day=0.0
    !status=1
    !return
    call GO_Timer_End( itim_write_day, status )
    IF_NOTOK_RETURN(status=1)


    call goLabel()
    status = 0
  end subroutine output_aerchemmip_write_daily


  subroutine output_aerchemmip_write_hourly(region,newhour,status)
    use MeteoData,      only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat   !
    use GO,             only : TDate, NewDate!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    use partools    , only : isRoot,myid
    implicit none
    logical,intent(in) ::newhour
    integer,intent(out)::status
    integer:: j_var
    integer::region
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_write_hourly'
    real                          :: rangehr,range6hr ! hour in days for bounds in NC-file

    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    ! root writes netcdf arrays
    real, allocatable :: arr3d(:,:,:) , arr2d(:,:)
    call goLabel(rname)

    call GO_Timer_Start( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
      write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!'
      call goErr
      status=1
      return
    end if

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

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

    ! allocate 3D and 4D global arrays for gathering data
    if (isRoot) then
      allocate( arr3d(imr,jmr,lmr)       )
      allocate( arr2d(imr,jmr)       )
    else
      ! other than root need the variable, but no space
      allocate( arr3d(1,1,1)   )
      allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 

    call date2tau( time_reftime6, itauref )
    ! call date2tau( idater, itaucur )

    rangehr=1.0/24.0!(3600)/86400.0
    do j_var =1,n_vars

      if (allvars(j_var)%table_id/='AERhr' .and. allvars(j_var)%table_id/='AER6hr' )then 
        cycle
      end if
      if (allvars(j_var)%dims==2)then
        if ( trim(allvars(j_var)%table_id)=='AERhr') then
          reftime = (itau - itauref) / 86400. - (1./24)

          allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_hr)
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)

          call gather( dgrid(1), allvars(j_var)%data2D , arr3d(:,:,1), 0, status)
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)               

          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d(:,:,1), status, start=(/1,1,timeidx_hr/), count=(/imr,jmr,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)


          ! write the date for timestep
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+(0.5/24.0)/) , status, start=(/timeidx_hr/), count=(/1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+rangehr/) , status, start=(/1,timeidx_hr/), count=(/2,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          ! Zero out the accumulated and written data for the next interval
          allvars(j_var)%data2D(i1:i2,j1:j2)=0.0
        end if
      else

        if ( trim(allvars(j_var)%table_id)=='AERhr') then
          reftime = (itau - itauref) / 86400. - (1./24)
          allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_hr)
          call GO_Timer_Start( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
          call GO_Timer_End( itim_write_gather, status )
          IF_NOTOK_RETURN(status=1)
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d, status, start=(/1,1,1,timeidx_hr/), count=(/imr,jmr,lmr,1/) )

          ! write the date for timestep
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime+(0.5/24.0)/) , status, start=(/timeidx_hr/), count=(/1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
          if (isroot) call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%bounds_varid,(/ reftime,reftime+(1./24.)/) , status, start=(/1,timeidx_hr/), count=(/2,1/) )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

          ! Zero out the accumulated and written data for the next interval
          allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0
        end if
      end if

      !end if

    end do
    deallocate(arr3d, arr2d)
    ! changed index to next hour
    timeidx_hr= timeidx_hr + 1
    ! accumulated hours to zero, actually this will always be 1h
    accumulation_hr=0.0


    !status=1
    !return
    call GO_Timer_End( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)


    call goLabel()
    status = 0
  end subroutine output_aerchemmip_write_hourly
  
  subroutine output_aerchemmip_write_6hourly(region,newhour,status)
    use MeteoData,      only : temper_dat,sst_dat,albedo_dat,ci_dat,lsp_dat,cp_dat,ssr_dat,str_dat,blh_dat,gph_dat,lwc_dat,iwc_dat,cco_dat,cc_dat,humid_dat, m_dat,phlb_dat,sp_dat   !
    use GO,             only : TDate, NewDate!
    use datetime,       only : tau2date,date2tau,julday   !
    use dims,           only : itau,iyear0          !current time
    use tm5_distgrid, only : dgrid,    Get_DistGrid ,gather
    use partools    , only : isRoot,myid
    use ebischeme,      only : AC_diag_prod,iprod_soa2dhour
    implicit none
    logical,intent(in) ::newhour
    integer,intent(out)::status
    integer::region
    integer:: j_var
    integer:: imr,jmr,i,i1,i2,j1,j2,lmr
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_write'

    ! reference time:
    integer, parameter            ::  time_reftime6(6) = (/1750,01,01,00,00,00/)
    integer(kind=8)                       ::  itauref  ! reftime in seconds
    real                          ::  reftime  ! seconds from reference time
    ! root writes netcdf arrays
    real, allocatable :: arr3d(:,:,:) , arr2d(:,:)
    call goLabel(rname)

    call GO_Timer_Start( itim_write_hour, status )
    IF_NOTOK_RETURN(status=1)

    if (region > 1) then
      write(gol,*) 'output_aerhemmip_write: region >1, only available in global mode!'
      call goErr
      status=1
      return
    end if

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

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

    ! allocate 3D and 4D global arrays for gathering data
    if (isRoot) then
      allocate( arr3d(imr,jmr,lmr)       )
      allocate( arr2d(imr,jmr)       )
    else
      ! other than root need the variable, but no space
      allocate( arr3d(1,1,1)   )
      allocate( arr2d(1,1)   )
    endif
    arr2d(:,:)=0.0
    arr3d(:,:,:)=0.0

    ! define the reference time in seconds (itauref)
    ! (for now in days since 1850-01-01 00:00, local variable) 

    call date2tau( time_reftime6, itauref )
    ! call date2tau( idater, itaucur )
    reftime = (itau - itauref) / 86400. 

    do j_var =1,n_vars

      if ( allvars(j_var)%table_id/='AER6hr' )then 
        cycle
      end if
      if (allvars(j_var)%dims==2)then
        !allvars(j_var)%data2D(i1:i2,j1:j2)=allvars(j_var)%data2D(i1:i2,j1:j2)/real(accumulation_6hr)
        call GO_Timer_Start( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)
        call gather( dgrid(1), allvars(j_var)%data2D , arr2d(:,:), 0, status)
        IF_NOTOK_RETURN(status=1)
        call GO_Timer_End( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr2d, status, start=(/1,1,timeidx_6hr/), count=(/imr,jmr,1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        ! write the date for timestep
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime/) , status, start=(/timeidx_6hr/), count=(/1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        ! Zero out the accumulated and written data for the next interval
        allvars(j_var)%data2D(i1:i2,j1:j2)=0.0

      else

        !allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)/real(accumulation_6hr)
        call GO_Timer_Start( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)
        call gather( dgrid(1), allvars(j_var)%data3D , arr3d(:,:,:), 0, status)
        call GO_Timer_End( itim_write_gather, status )
        IF_NOTOK_RETURN(status=1)
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%varid,arr3d, status, start=(/1,1,1,timeidx_6hr/), count=(/imr,jmr,lmr,1/) )

        ! write the date for timestep
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,allvars(j_var)%time_varid,(/ reftime/) , status, start=(/timeidx_6hr/), count=(/1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        !start=(/i1,j1,1,timeidx_mon/), count=(/imr,jmr,lmr,1/) )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        ! Zero out the accumulated and written data for the next interval
        allvars(j_var)%data3D(i1:i2,j1:j2,1:lmr)=0.0

      end if
    end do
    deallocate(arr3d,arr2d)
    ! changed index to next 6hour
    timeidx_6hr= timeidx_6hr + 1
    ! exception for one 6hr field, no need to do another subroutine for it
    accumulation_6hr=0.0

    ! zero out the chemical production (SOA for hourly output)
    ! now used fro 6h output
    !AC_diag_prod(region)%prod(i1:i2,j1:j2,:,iprod_soa2dhour)=0.0

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

    call goLabel()
    status = 0
  end subroutine output_aerchemmip_write_6hourly

  subroutine accumulate_data(dhour,l_do_ec550aer_only,status)
    use GO         , only : TDate, NewDate, rTotal, operator(-)
    use Grid       , only : FPressure,HPressure
    use binas,       only : rgas, rol,xmair,Dobs,Avog
    USE toolbox,       only : ltropo_ifs,     lvlpress
    !use datetime,    only : tau2date
    use MeteoData,      only : temper_dat, sst_dat, albedo_dat, ci_dat, lsp_dat, cp_dat, ssr_dat, str_dat, blh_dat, &
         gph_dat, lwc_dat, iwc_dat, cco_dat, cc_dat, humid_dat, m_dat, phlb_dat, sp_dat, pu_dat, pv_dat,pw_dat
    use photolysis_data,only:phot_dat   !
    use global_data, only : mass_dat, region_dat,conv_dat
    use tracer_data, only : chem_dat
    use dims,        only : lm,sec_month
    use chem_param,  only : iso4acs, iso4cos, iduacs, iso4ais, ibccos, ibcaii, xmair, xmo3,nmhc,xmcb5,ncb5, o3p, o3l,ino3_a,xmn,ra
    !use chem_param,  only : iso4nus, isscos, ino3_a, issacs,  iducos, iduaci, nmod
    !use chem_param,  only : iducoi, ibcacs, ipomcos, ipomaii, ibcais, ipomacs, ipomais
    !use chem_param,  only : imsa, inh4
    !use chem_param,  only : ntrace,names,mode_end_so4
    !use mo_aero_m7,  only : nmod!,nsol
    !use optics,      only : optics_aop_get
    use m7_data,     only : h2o_mode, rw_mode, rwd_mode
    use sedimentation, only : nsed,ised ,sindex
    use sedimentation, only : deposition => buddep_m7_dat   !(i,j,lev,nsed)
    use wet_deposition,only : wetdep=>buddep_dat            !(i,j,lev,ntrace)
    use emission_data,  only     : budemi_dat
    use ebischeme,      only : buddry_dat => buddep_dat
    use chem_param,  only : xmso2, xmso4, xmdms, xmpom, xmbc, xmdust, xmnacl,xmo3,xmnh4,ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr,ra
    use calc_pm,     only : PMx_Integrate_3d,PMx_integrate_0d
    use emission_nox,  only : eminox_lightning
    use ebischeme,   only : diag_prod,AC_diag_prod,AC_O3_lp,AC_loss,iloss_ch4,iloss_co,iprod_gasso4,iprod_aqso4,iprod_soa3dmon,iprod_soa2dhour
    use partools    , only : isRoot,myid
    use dims,    only: gtor, dx, dy, ybeg, xref, yref,ndyn
    use binas,   only: ae
    use emission_nox, only: flashrate_lightning
    implicit none

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

    !    integer                 :: indices(7)
    integer                 :: itrac,gasind
    integer                 :: i_sed
    integer                 :: i_emi
    integer                 :: i, j, k, imr, jmr, lmr, lwl, dtime,index,imode,m
    integer                 :: i1, i2, j1, j2
    integer                 :: ie, je ! indices for subdomain extended with halo cells
    integer,intent(in)      :: dhour
    integer                 :: status
    integer                 :: j_var,region,i_var,i_wdep,sedindex,icomp

    real                    :: dens
    real                    :: vol

    real                    :: tempbud,xmcomp,temp
    real, dimension(:,:,:,:), pointer :: tracers     ! transported tracers
    real, dimension(:,:,:,:), pointer :: tracers_c   ! non-transported
    real, dimension(:),       pointer :: dxyp
    integer, dimension(n_indices) :: indices
    real::xmgas
    real, dimension(:,:,:), allocatable :: pres3d
    real, dimension(:,:,:), allocatable :: temp_pm
    real, dimension(:,:,:), allocatable :: pres3dh
    real :: pm_sizelimit
    integer::emi_index,wet_index,dry_index

    ! tropopause calculations
    real, dimension(:), allocatable :: gphx, tx
    real, dimension(:,:,:), pointer :: gph   ! height (incl. oro)
    real, dimension(:,:,:), pointer :: t     ! temperature (K)
    integer :: itropo
    real::xres,yres,dxx,dyy,uwind,vwind,lat,wwind,dz,meanwind
    integer, dimension(4) :: modes_dust=(/3,4,6,7/)
    !
    !EC550AER
    logical,intent(in) :: l_do_ec550aer_only

    call goLabel(rname)

    call GO_Timer_Start( itim_accu, status )
    IF_NOTOK_RETURN(status=1)
    region=1
    ! 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

    ! for tropopause variables
    allocate(gphx(0:lm(region)))   ! note now from 0-->lm
    allocate(tx(lm(region)))
    t           => temper_dat (region)%data
    gph         =>    gph_dat (region)%data

    allocate( temp_pm(i1:i2,j1:j2,lmr) )
    temp_pm=0.0
    allocate( pres3d(i1:i2,j1:j2,lmr) )
    allocate( pres3dh(i1:i2,j1:j2,0:lmr) )

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



    accumulation_6hr=0.0!accumulation_6hr+dtime


    ! Gridbox area
    dxyp   => region_dat(region)%dxyp

    ! mass_dat and chem_data keep data in kg/gridbox
    !
    tracers => mass_dat(region)%rm
    tracers_c => chem_dat(region)%rm

    if (.not. l_do_ec550aer_only) then
      dtime = dhour*3600
      accumulation_mon=accumulation_mon+dtime
      accumulation_hr=accumulation_hr+dtime
      accumulation_day=accumulation_day+dtime
      do j_var = 1, n_vars
        indices(:)=allvars(j_var)%indices(:)
        !Here we use the tm5-indices to collect data for output
        !

        if (trim(allvars(j_var)%class)=='meteo2D') then
          do j = j1,j2
            do i=i1,i2
              tx(:)=t(i,j,:)
              gphx(0:lm(region))=gph(i,j,1:lm(region)+1)   !note the bounds
              !ibase  = lubottom(i,j)
              !itop   = lutop(i,j)
              itropo = ltropo_ifs(region,i,j,tx,lm(region))
              ! 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)
              if (trim(allvars(j_var)%vname)=='ps')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)
              else if  (trim(allvars(j_var)%vname)=='ptp')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+pres3d(i,j,itropo)*dtime
              else if (trim(allvars(j_var)%vname)=='tropoz')then
                ! conversion (in order of execution):
                !kg->g (1e3), g->molec (xmo3), m2->cm2(1e-4) , molec-> molec/m2(dxyp) , molec/cm2->dobson (DOBS)
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*sum(tracers(i,j,1:itropo,io3)*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs)
              else if (trim(allvars(j_var)%vname)=='toz')then
                ! conversion (in order of execution):
                !kg->g (1e3), g->molec (xmo3), m2->cm2(1e-4) , molec-> molec/m2(dxyp) , molec/cm2->dobson (DOBS)
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*sum(tracers(i,j,:,io3)*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs)
              else if  (trim(allvars(j_var)%vname)=='tos')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)!+dtime*sst_dat(1)%data(i,j,1)
              else if  (trim(allvars(j_var)%vname)=='sic')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*ci_dat(1)%data(i,j,1)
              else if  (trim(allvars(j_var)%vname)=='tatp')then
                ! t at the gridpoint center ->mean at the interface value
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*t(i,j,itropo)
              else if  (trim(allvars(j_var)%vname)=='ztp')then
                !gph at interface
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(gph(i,j,itropo+1)+gph(i,j,itropo))/2
              else if  (trim(allvars(j_var)%vname)=='bldep')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*conv_dat(1)%blh(i,j)
              else if  (trim(allvars(j_var)%vname)=='pr')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(lsp_dat(1)%data(i,j,1)+ cp_dat(1)%data(i,j,1))
              else if  (trim(allvars(j_var)%vname)=='albsrfc')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*albedo_dat(1)%data(i,j,1)
              end if
            end do
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='meteo3D') then
          do j = j1,j2
            do i=i1,i2
              if (trim(allvars(j_var)%vname)=='phalf')then
                do k=1,lmr+1
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*pres3dh(i,j,k-1)
                end do
              else
                do k =1,lmr
                  if(trim(allvars(j_var)%vname)=='ta')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*temper_dat(1)%data(i,j,k)
                  else if  (trim(allvars(j_var)%vname)=='emilnox')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime* eminox_lightning(region)%d3(i,j,k) / dxyp(j)
                  else if  (trim(allvars(j_var)%vname)=='jno2')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*phot_dat(region)%jno2(i,j,k)
                  else if  (trim(allvars(j_var)%vname)=='airmass')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime* m_dat(region)%data(i,j,k) / dxyp(j)
                  else if  (trim(allvars(j_var)%vname)=='pfull')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*pres3d(i,j,k)
                    !                      else if  (trim(allvars(j_var)%vname)=='phalf')then
                    !                         allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*pres3dh(i,j,k)
                  else if  (trim(allvars(j_var)%vname)=='hus')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*humid_dat(1)%data(i,j,k)
                  else if  (trim(allvars(j_var)%vname)=='zg')then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*gph_dat(1)%data(i,j,k)
                  else if  (trim(allvars(j_var)%vname)=='ua')then
                    ! distance of single gridbox
                    yres = dy/yref(1)
                    xres = dx/xref(1)
                    lat = ybeg(1) + 0.5 * yres + yres * (j-1)
                    dxx = ae * xres * gtor * cos(lat*gtor)
                    uwind=dxx*(pu_dat(1)%data(i,j,k) + pu_dat(1)%data(i-1,j,k))*0.5 / m_dat(1)%data(i,j,k)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*uwind
                  else if  (trim(allvars(j_var)%vname)=='va')then
                    yres = dy/yref(1)
                    dyy = ae * yres * gtor
                    vwind= dyy   *(pv_dat(1)%data(i,j,k) + pv_dat(1)%data(i,j+1,k))*0.5 / m_dat(1)%data(i,j,k)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*vwind
                  else if  (trim(allvars(j_var)%vname)=='wa')then
                    dz = gph_dat(1)%data(i,j,k+1) - gph_dat(1)%data(i,j,k)
                    wwind= dz*(pw_dat(1)%data1(i,j,k-1) + pw_dat(1)%data1(i,j,k))*0.5 / m_dat(1)%data(i,j,k)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*wwind
                  end if
                end do
              end if
            end do
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='emon') then
          !Sedimented indices in deposition%sed
          !do i_sed=1,nsed
          !gridpoints
          do j = j1,j2
            do i=i1,i2
              select case ( trim(allvars(j_var)%vname))
              case('loaddust')
                do m=1,n_indices
                  index=allvars(j_var)%indices(m)
                  if (index==0) cycle
                  do k=1,lmr
                    allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,k,index))/ dxyp(j)
                  end do
                end do
              case('concdust')
                do m=1,n_indices
                  index=allvars(j_var)%indices(m)
                  if (index==0) cycle
                  do k=1,lmr
                    vol = (gph_dat(region)%data(i,j,k+1)-gph_dat(region)%data(i,j,k)) * dxyp(j)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*(tracers(i,j,k,index))/ vol
                  end do
                end do
              case('seddustCI')
                sedindex=sindex(iducoi)
                xmcomp=xmdust
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+ &
                     deposition(region)%sed(i,j,1,sedindex) / dxyp(j) * 1.E-03*xmcomp*-100000
              case('depdust')
                tempbud=0.0
                xmcomp=xmdust
                do m=1,n_indices
                  !exit if indices ==0
                  !since indices after 0 will all be also 0
                  if (indices(m)>0) then
                    if (allvars(j_var)%indices(m)>69)then
                      sedindex=-1
                    else
                      sedindex=sindex(allvars(j_var)%indices(m))
                    end if
                    if (sedindex>0)then
                      tempbud=tempbud+( &
                           sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp
                    end if
                  end if
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud)
                end do
              case('md')
                do k=1,lmr
                  ! tendency of atm_mass content of dust dry aerosol particles due to emission
                  ! coarse mode only 3d hmmm
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)-1!dtime*rw_mode(1,7)%d3(i,j,k)
                end do
              case('flashrate')
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*flashrate_lightning(1)%d23(i,j)/(dxyp(j)*1e-4) !km-2s-1  (flashrate_lightning in [s-1]
              case('conccn')
                do k=1,lmr
                  temp=0.0
                  dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas)
                  do  m=1,n_indices
                    index=allvars(j_var)%indices(m)
                    if (index==0) cycle
                    temp=tracers(i,j,k,index)*dens
                  end do
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*temp
                end do
              case('sconcss')
                temp=0.0
                k=1 !Surface
                dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas)
                do  m=1,n_indices
                  index=allvars(j_var)%indices(m)
                  if (index==0) cycle
                  temp=tracers(i,j,k,index)*dens
                end do
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*temp
              case DEFAULT !('default')
                write(gol,'("emon accumulate missing: ",a)') trim(allvars(j_var)%vname); call goPr
              end select
            end do
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='dry') then
          !Sedimented indices in deposition%sed
          !do i_sed=1,nsed
          !gridpoints
          do j = j1,j2
            do i=i1,i2
              select case ( trim(allvars(j_var)%vname))
              case('dryoa')
                dry_index=1
                xmcomp=xmpom
              case('drybc')
                dry_index=2
                xmcomp=xmbc
              case ('dryso2')
                dry_index=3
                xmcomp=xmso2
              case('dryso4')
                dry_index=4
                xmcomp=xmso4
              case('drydust')
                dry_index=5
                xmcomp=xmdust
              case('drydms')
                dry_index=6
                xmcomp=xmdms
              case('dryss')
                dry_index=7
                xmcomp=xmnacl
              case('drysoa')
                dry_index=8
                xmcomp=xmpom!soa
              case('drynh3')
                dry_index=9
                xmcomp=xmnh3
              case('drynh4')
                dry_index=10
                xmcomp=xmnh4
              case('drynoy')
                dry_index=11
                xmcomp=1
              case('dryo3')
                dry_index=12
                xmcomp=xmo3
              case('dryno3')
                xmcomp=xmno3
              case('dryno2')
                xmcomp=xmno2
              case('dryno')
                xmcomp=xmno
              case DEFAULT !('default')
                write(gol,'("dry xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr
              end select
              tempbud=0.0
              !if (xx='dryoa') then
              !end if
              do m=1,n_indices
                !exit if indices ==0
                !since indices after 0 will all be also 0
                if (indices(m)<=0) then
                  cycle
                else
                  if (indices(m)>69) then
                    sedindex=-1
                  else
                    sedindex=sindex(allvars(j_var)%indices(m))
                  end if
                  !only M7 aerosol tracers (and msa/nh4/no3) are sedimented and deposited
                  if ( trim(allvars(j_var)%vname)=='drynoy')then
                    select case(allvars(j_var)%indices(m))
                    case(ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr)
                      xmcomp=xmn
                    case(in2o5)
                      xmcomp=xmn*2.0
                    end select
                  end if
                  if (sedindex>0)then
                    tempbud=tempbud+(buddry_dat(region)%dry(i,j,indices(m)) +&
                         sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                    !for others only deposition applies
                  else
                    tempbud=tempbud+buddry_dat(region)%dry(i,j,allvars(j_var)%indices(m)) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                  end if
                end if
              end do
              !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-drydepos(region)%f2dslast(i,j,dry_index))
              !drydepos(region)%f2dslast(i,j,dry_index)=tempbud
              allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
              allvars(j_var)%budgetdata(i,j,1)=tempbud
            end do
          end do

          !end if
        else if (trim(allvars(j_var)%class)=='wet') then
          !gridpoints
          do j = j1,j2
            do i=i1,i2
              select case ( trim(allvars(j_var)%vname))
              case('wetoa')
                wet_index=1
                xmcomp=xmpom
              case('wetbc')
                wet_index=2
                xmcomp=xmbc
              case ('wetso2')
                wet_index=3
                xmcomp=xmso2
              case('wetso4')
                wet_index=4
                xmcomp=xmso4
              case('wetdust')
                wet_index=5
                xmcomp=xmdust
              case('wetnoy')
                !wet_index=6 !DMS
                wet_index=13
                xmcomp=1
              case('wetss')
                wet_index=7
                xmcomp=xmnacl
              case('wetsoa')
                wet_index=8
                xmcomp=xmpom!soa
              case('wetnh3')
                wet_index=9
                xmcomp=xmnh3
              case('wetnh4')
                wet_index=10
                xmcomp=xmnh4
              case DEFAULT !('default')
                write(gol,'("wetdep xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr
              end select
              tempbud=0.0
              if ( trim(allvars(j_var)%vname)=='wetnoy')then
                do m=1,n_indices
                  select case(allvars(j_var)%indices(m))
                  case(ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr)
                    xmcomp=xmn
                  case(in2o5)
                    xmcomp=xmn*2.0
                  end select
                  if(allvars(j_var)%indices(m)>0 .and. allvars(j_var)%indices(m)<=ntracet)then
                    tempbud=tempbud+(sum(wetdep(region)%lsp(i,j,:,indices(m))) +   sum(wetdep(region)%cp(i,j,:,indices(m)))) / dxyp(j) * 1.E-03*xmcomp
                  end if
                end do
                !write(gol,'("wet: ",a,i3,2E30.7e5)') trim(allvars(j_var)%vname),indices(m),sum(wetdep(region)%lsp(i,j,:,indices(m)))/ dxyp(j) * 1.E-03*xmcomp, sum(wetdep(region)%cp(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp; call goPr
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) !kg m-2s-1 (s-1 at the time of writing)
                allvars(j_var)%budgetdata(i,j,1)=tempbud
              else
                do m=1,n_indices
                  !exit if indices ==0
                  !since indices after 0 will all be also 0, and wetdep not defined for non-transported vars
                  if (indices(m)<=0.or.indices(m)>ntracet) then
                    cycle
                  else
                    if ( trim(allvars(j_var)%vname)=='wetnoy')then
                      select case(indices(m))
                      case(ino,ino2,ino3,ihno3,ino3_a,ihno4,ihono,ich3o2no2,ipan,iorgntr)
                        xmcomp=xmn
                      case(in2o5)
                        xmcomp=xmn*2.0
                      end select                      ! go through gridpoints
                    end if
                    tempbud=tempbud+(sum(wetdep(region)%lsp(i,j,:,indices(m))) +   sum(wetdep(region)%cp(i,j,:,indices(m)))) / dxyp(j) * 1.E-03*xmcomp
                  end if
                end do
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1)) !kg m-2s-1 (s-1 at the time of writing)
                allvars(j_var)%budgetdata(i,j,1)=tempbud
              end if
            end do
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='che') then
          do m=1,n_indices
            !exit if indices ==0
            !since indices after 0 will all be also 0
            if (indices(m)<=0) then
              cycle
            else
              ! go through gridpoints
              do j = j1,j2
                do i=i1,i2
                  do k=1,lmr
                    ! add emitted mass from different size ranges (i_emi)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime! / dxyp(j)) * 1.E-03
                  end do
                end do
              end do
            end if
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='chep') then
          ! go through gridpoints
          do j = j1,j2
            do i=i1,i2
              do k=1,lmr
                ! density for conversion of aerosol mass densities ( ---> 1/m3 )
                ! (conversion factor 1.E-03 is for  g --> kg)
                !volume(i1:ie,j,1:lmr) = (gph_dat(region)%data(i1:ie,j,2:lmr+1)-gph_dat(region)%data(i1:ie,j,1:lmr)) * dxyp(j)
                vol = (gph_dat(region)%data(i,j,k+1)-gph_dat(region)%data(i,j,k)) * dxyp(j)
                ! dtime not needed as time is already included in calculation of prod-field
                ! gas-phase so4 production
                if (trim(allvars(j_var)%vname)=='chegpso4')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_diag_prod(region)%prod(i,j,k,iprod_gasso4)/ dxyp(j)
                  !liquid so4 production
                else if (trim(allvars(j_var)%vname)=='cheaqpso4')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_diag_prod(region)%prod(i,j,k,iprod_aqso4)/ dxyp(j)
                else if (trim(allvars(j_var)%vname)=='lossch4')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_loss(region)%prod(i,j,k,iloss_ch4)/vol
                else if (trim(allvars(j_var)%vname)=='lossco')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_loss(region)%prod(i,j,k,iloss_co)/vol
                else if (trim(allvars(j_var)%vname)=='o3loss')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_O3_lp(region)%prod(i,j,k,1) / vol
                else if (trim(allvars(j_var)%vname)=='o3prod')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_O3_lp(region)%prod(i,j,k,2) /vol
                end if
              end do
              !if (trim(allvars(j_var)%vname)=='chepsoa')then
              !   !All soa (svoc+elvoc)
              !   allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa3dmon))/ dxyp(j)
              !end if
            end do
          end do
          if (trim(allvars(j_var)%vname)=='chepsoa')then
            do j =j1,j2
              do i =i1,i2
                !All soa (svoc+elvoc)
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa3dmon))/ dxyp(j)
              end do
            end do
          end if

          !end if
        else if (trim(allvars(j_var)%class)=='emi') then
          !Data can be found in emis3d in emission_pom.f90, only when declaring emissions
          ! which requires to write from there.
          ! or in pom_emis_2/3d(region,sectors), available all the time.
          !Could also use mode based index holder: nmode_end_XXX(nmod), where XXX=[so4,bc,pom,ss,dust]
          !Sedimented indices in deposition%sed
          do j = j1,j2
            do i=i1,i2
              ! add emitted mass from different size ranges (i_emi)
              select case ( trim(allvars(j_var)%vname))
              case('emioa')
                emi_index=1
                xmcomp=xmpom
              case('emibc')
                emi_index=2
                xmcomp=xmbc
              case ('emiso2')
                emi_index=3
                xmcomp=xmso2
              case('emiso4')
                emi_index=4
                xmcomp=xmso4
              case('emidust')
                emi_index=5
                xmcomp=xmdust
              case('emidms')
                emi_index=6
                xmcomp=xmdms
              case('emiss')
                emi_index=7
                xmcomp=xmnacl
              case('emiisop')
                emi_index=8
                xmcomp=xmisop
              case('emivoc')
                emi_index=9
                xmcomp=1!voc
              case('eminh3')
                emi_index=10
                xmcomp=xmnh3
              case('eminh4')
                emi_index=11
                xmcomp=xmnh4
              case('emico')
                emi_index=12
                xmcomp=xmco
              case('emibvoc')
                emi_index=13
                xmcomp=xmterp
              case('eminox')
                !emi_index=13
                xmcomp=xmn
              case('emipoa')
                emi_index=1
                xmcomp=xmpom
              case DEFAULT
                write(gol,'("emi xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr
              end select
              tempbud=0.0
              if (trim(allvars(j_var)%vname)=='emivoc')then
                do m=1,ncb5
                  select case(nmhc(m))
                  case(ipar,ich2o,ich3oh,ihcooh)
                    xmcomp=xmc
                  case(ieth,iole,iald2,imcooh,ic2h6,iethoh)
                    xmcomp=xmc*2.0
                  case(imgly,ic3h8,ic3h6,iacet)
                    xmcomp=xmc*3.0
                  end select
                  tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,nmhc(m))) / dxyp(j) * 1.E-03*xmcomp!kg(C)m-2s-1
                end do
                do m=1,2
                  if (m==1)then
                    gasind=iterp
                    xmcomp=xmc*5*2
                  else
                    gasind=iisop
                    xmcomp=xmc*5.0
                  end if
                  tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp!kg(C)m-2s-1
                end do
                !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-emission(region)%f2dslast(i,j,emi_index))
                !emission(region)%f2dslast(i,j,emi_index)=tempbud
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                allvars(j_var)%budgetdata(i,j,1)=tempbud
              else if (trim(allvars(j_var)%vname)=='emibvoc')then
                tempbud=0.0
                do m=1,2!ncb5
                  if (m==1)then
                    gasind=iterp
                    xmcomp=2*5*xmc!xmterp=2*cxmisop and as C: 2*5*xmc
                  else
                    gasind=iisop
                    xmcomp=5*xmc!xmisop as mass oc C
                  end if
                  tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp
                end do
                !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-emission(region)%f2dslast(i,j,emi_index))
                !emission(region)%f2dslast(i,j,emi_index)=tempbud
                allvars(j_var)%data2D(i,j)= 0.0!allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                allvars(j_var)%budgetdata(i,j,1)=tempbud
              else if (trim(allvars(j_var)%vname)=='eminox') then
                tempbud=0.0
                indices=allvars(j_var)%indices
                do m=1,n_indices
                  if (indices(m)<=0)then
                    cycle
                  else
                    xmcomp=xmn
                    !write(gol,'("eminox: ",voca,i3,E20.7e1)')trim(allvars(j_var)%vname), indices(m),sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp; call goPr
                    tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp !kg(N) m-2s-1  (s-1 at the time of writing)
                  end if
                end do
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                allvars(j_var)%budgetdata(i,j,1)=tempbud
              else
                tempbud=0.0
                ! add soa for emioa (emission + production)
                do m=1,n_indices
                  !exit if indices ==0
                  !since indices after 0 will all be also 0
                  if (indices(m)<=0) then
                    cycle
                  else
                    ! go through gridpoints
                    tempbud = tempbud + sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp
                    allvars(j_var)%data2D(i,j) = allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  end if
                end do
              end if
              ! add soa production to oa emission after poa is calculated.
              if  (allvars(j_var)%vname == 'emioa' .and. trim(allvars(j_var)%freq)=='mon') then
                allvars(j_var)%data2D(i,j) = allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,3))/ dxyp(j)
              end if
            end do
          end do

          !end if
          ! CRESCENDO extensions !
        else if (trim(allvars(j_var)%class)=='crescendo')then
          if (trim(allvars(j_var)%vname)=='od5503ddust')then
            cycle
          else if (trim(allvars(j_var)%vname)=='od550dust')then
            cycle
          else if (trim(allvars(j_var)%vname)=='od550aer')then
            cycle
          else if (trim(allvars(j_var)%vname)=='ec550aer')then
            cycle
          else if (trim(allvars(j_var)%vname)=='od440aer')then
            cycle
          else if (trim(allvars(j_var)%vname)=='od440dust')then
            cycle
          else if (trim(allvars(j_var)%vname)=='od870aer')then
            cycle
          else if (trim(allvars(j_var)%vname)=='dmsos')then
            cycle
          END if
          index=allvars(j_var)%indices(1)
          indices=allvars(j_var)%indices(:)
          if (allvars(j_var)%dims==3)then
            do k=1,lmr
              do j=j1,j2
                do i=i1,i2
                  if (trim(allvars(j_var)%vname)=='ndtot')then
                    dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas)
                    do m=1,7
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*(tracers(i,j,k,indices(m)))*dens
                    end do
                  else if  (trim(allvars(j_var)%vname)=='emilnox')then
                    ! eminoc_lightning is /month so dtime/sec_month gives us per this step
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime/sec_month* eminox_lightning(region)%d3(i,j,k) / dxyp(j)
                  else if (trim(allvars(j_var)%vname)=='mmraerh2o_1')then
                    allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,1)%d3(i,j,k))/m_dat(region)%data(i,j,k)
                  else if (trim(allvars(j_var)%vname)=='mmraerh2o_2')then
                    allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,2)%d3(i,j,k))/m_dat(region)%data(i,j,k)
                  else if (trim(allvars(j_var)%vname)=='mmraerh2o_3')then
                    allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,3)%d3(i,j,k))/m_dat(region)%data(i,j,k)
                  else if (trim(allvars(j_var)%vname)=='mmraerh2o_4')then
                    allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,4)%d3(i,j,k))/m_dat(region)%data(i,j,k)
                  else if (trim(allvars(j_var)%vname)=='chepsoa3d')then
                    !All soa (svoc+elvoc)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+AC_diag_prod(region)%prod(i,j,k,iprod_soa3dmon)/ dxyp(j)
                    !number concentrations (ndtot, sfndtot and sfnd100 handled separately)
                  else if (trim(allvars(j_var)%unit)==units(im3unit).and. .not. (trim(allvars(j_var)%vname)=='sfndtot' .or. trim(allvars(j_var)%vname)=='sfnd100' .or. trim(allvars(j_var)%vname)=='ndtot' )   ) then
                    dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas)
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*dens
                  else if (trim(allvars(j_var)%unit)==units(im3unit)) then
                    !nd? output processed here
                    dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas)
                    !1e-6 for m3->cm-3
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*dens
                  else if (index>0 .and. index<=ntracet) then
                    if  (trim(allvars(j_var)%unit)==units(ivmrunit)) then
                      xmcomp=allvars(j_var)%xmgas
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*xmair/xmcomp/m_dat(region)%data(i,j,k) !kmole kmole-1
                    else if  (trim(allvars(j_var)%unit)==units(immrunit)) then
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1
                    end if
                  else if (index>ntracet)then ! than non transported (tracers_c=> chem_dat(region)%rm)
                    if  (trim(allvars(j_var)%unit)==units(ivmrunit)) then
                      xmcomp=allvars(j_var)%xmgas
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)*xmair/xmcomp/m_dat(region)%data(i,j,k)
                    else if  (trim(allvars(j_var)%unit)==units(immrunit)) then
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1
                    else
                      write(gol,'("user_output_aerchemmip: no case for output ",a)') trim(allvars(j_var)%vname) ; call goErr
                      TRACEBACK
                      status=1; return
                    end if
                  else
                    write(gol,*) 'user_output_aerchemmip: missing accumulation for CRESCENDO 2D variable ',allvars(j_var)%vname ; call goErr
                    TRACEBACK
                    status=1; return
                  end if
                end do
              end do
            end do
          else !2D vars
            do j=j1,j2
              do i=i1,i2
                tempbud=0.0
                if (trim(allvars(j_var)%vname)=='sfmmrss')then
                  tempbud=0.0
                  do m=1,size(allvars(j_var)%indices)
                    index=allvars(j_var)%indices(m)
                    if (index>0)then
                      tempbud=tempbud+tracers(i,j,1,index)/m_dat(region)%data(i,j,1)!kg kg-1
                    end if
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den
                  allvars(j_var)%data2D(i,j)= tempbud
                else if (trim(allvars(j_var)%vname)=='co2mass')then
                  !NOTAVAILABLE
                  do k=1,lmr
                    !mass of CO2, sum over levels
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers_c(i,j,k,iaco2))!/m_dat(region)%data(i,j,k)
                  end do
                else if (trim(allvars(j_var)%vname)=='chepsoa2d')then
                  !All soa (svoc+elvoc)
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa2dhour))/ dxyp(j)
                  allvars(j_var)%data2D(i,j)= sum(AC_diag_prod(region)%prod(i,j,:,iprod_soa2dhour))/ dxyp(j)/dtime !
!!! FOLLOWING CRESCENDO VARS: Instantaneous, no accumulation
                else if (trim(allvars(j_var)%vname)=='sfmmrso4')then
!!! Instantaneous, no accumulation
                  tempbud=0.0
                  do m=1,size(allvars(j_var)%indices)
                    index=allvars(j_var)%indices(m)
                    if (index>0)then
                      !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den
                      tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1
                    end if
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den
                  allvars(j_var)%data2D(i,j)= tempbud
                else if (trim(allvars(j_var)%vname)=='sfmmrbc')then
!!! Instantaneous, no accumulation
                  tempbud=0.0
                  do m=1,size(allvars(j_var)%indices)
                    index=allvars(j_var)%indices(m)
                    if (index>0)then
                      tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1
                    end if
                  end do
                  allvars(j_var)%data2D(i,j)=tempbud
                else if (trim(allvars(j_var)%vname)=='sfmmrsoa')then
!!! Instantaneous, no accumulation
                  tempbud=0.0
                  do m=1,size(allvars(j_var)%indices)
                    index=allvars(j_var)%indices(m)
                    if (index>0)then
                      !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den
                      tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1
                    end if
                  end do
                  allvars(j_var)%data2D(i,j)= tempbud
                else if (trim(allvars(j_var)%vname)=='sfmmroa')then
!!! Instantaneous, no accumulation
                  tempbud=0.0
                  do m=1,size(allvars(j_var)%indices)
                    index=allvars(j_var)%indices(m)
                    if (index>0)then
                      !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den
                      tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1
                    end if
                  end do
                  allvars(j_var)%data2D(i,j)= tempbud
                else if (trim(allvars(j_var)%vname)=='sfisop')then
!!! Instantaneous, no accumulation
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,iisop)*xmair/xmisop/m_dat(region)%data(i,j,1)
                  allvars(j_var)%data2D(i,j)=tracers(i,j,1,iisop)*xmair/xmisop/m_dat(region)%data(i,j,1) !kilomole kilomole-1
                else if (trim(allvars(j_var)%vname)=='sfdms')then
!!! Instantaneous, no accumulation
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,idms)*xmair/xmdms/m_dat(region)%data(i,j,1)
                  allvars(j_var)%data2D(i,j)= tracers(i,j,1,idms)*xmair/xmdms/m_dat(region)%data(i,j,1) !kmole kmole-1
                else if (trim(allvars(j_var)%vname)=='sfmono')then
!!! Instantaneous, no accumulation
                  !tracers_c(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k)
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,iterp)*xmair/xmterp/m_dat(region)%data(i,j,1)
                  allvars(j_var)%data2D(i,j)= tracers(i,j,1,iterp)*xmair/xmterp/m_dat(region)%data(i,j,1) !kmole kmole-1
!!! no accumulation
                else if (trim(allvars(j_var)%vname)=='sfmmrdustabv1')then
!!! Instantaneous, no accumulation
                  allvars(j_var)%data2D(i,j)= (&
                       (1.0-mode_fraction(rw_mode(1,3)%d3(i,j,1),1000e-9,3))*tracers(i,j,1,iduacs)+&
                       (1.0-mode_fraction(rw_mode(1,4)%d3(i,j,1),1000e-9,4))*tracers(i,j,1,iducos)+&
                       (1.0-mode_fraction(rw_mode(1,6)%d3(i,j,1),1000e-9,6))*tracers(i,j,1,iduaci)+&
                       (1.0-mode_fraction(rw_mode(1,7)%d3(i,j,1),1000e-9,7))*tracers(i,j,1,iducoi))&
                       /m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfmmrdustabv10')then
!!! Instantaneous, no accumulation
                  allvars(j_var)%data2D(i,j)=(&
                       (1.0-mode_fraction(rw_mode(1,3)%d3(i,j,1),10000e-9,3))*tracers(i,j,1,iduacs)+&
                       (1.0-mode_fraction(rw_mode(1,4)%d3(i,j,1),10000e-9,4))*tracers(i,j,1,iducos)+&
                       (1.0-mode_fraction(rw_mode(1,6)%d3(i,j,1),10000e-9,6))*tracers(i,j,1,iduaci)+&
                       (1.0-mode_fraction(rw_mode(1,7)%d3(i,j,1),10000e-9,7))*tracers(i,j,1,iducoi))&
                       /m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfmmrdustbel1')then
!!! Instantaneous, no accumulation
                  allvars(j_var)%data2D(i,j)= (&
                       (mode_fraction(rw_mode(1,3)%d3(i,j,1),1000e-9,3))*tracers(i,j,1,iduacs)+&
                       (mode_fraction(rw_mode(1,4)%d3(i,j,1),1000e-9,4))*tracers(i,j,1,iducos)+&
                       (mode_fraction(rw_mode(1,6)%d3(i,j,1),1000e-9,6))*tracers(i,j,1,iduaci)+&
                       (mode_fraction(rw_mode(1,7)%d3(i,j,1),1000e-9,7))*tracers(i,j,1,iducoi))&
                       /m_dat(region)%data(i,j,1)!kg kg-1
                else if  (trim(allvars(j_var)%vname)=='uas')then
                  ! distance of single gridbox
                  yres = dy/yref(1)
                  xres = dx/xref(1)
                  lat = ybeg(1) + 0.5 * yres + yres * (j-1)
                  dxx = ae * xres * gtor * cos(lat*gtor)
                  uwind=dxx*(pu_dat(1)%data(i,j,1) + pu_dat(1)%data(i-1,j,1))*0.5 / m_dat(1)%data(i,j,1)
                  allvars(j_var)%data2D(i,j)= uwind
                else if  (trim(allvars(j_var)%vname)=='vas')then
                  yres = dy/yref(1)
                  dyy = ae * yres * gtor
                  vwind= dyy   *(pv_dat(1)%data(i,j,1) + pv_dat(1)%data(i,j+1,1))*0.5 / m_dat(1)%data(i,j,1)
                  allvars(j_var)%data2D(i,j)= +vwind
!!!6hr output ending here
                else if (trim(allvars(j_var)%vname)=='ps' .and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)=sp_dat(1)%data(i,j,1)!Pa
                else if (trim(allvars(j_var)%vname)=='sfno'.and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)=tracers_c(i,j,1,ino)*xmair/xmno/m_dat(region)%data(i,j,1)!mole mole-1
                else if (trim(allvars(j_var)%vname)=='sfnh3'.and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)= tracers(i,j,1,inh3)*xmair/xmnh3/m_dat(region)%data(i,j,1)!mole mole-1
                else if (trim(allvars(j_var)%vname)=='sfhno3'.and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)= tracers(i,j,1,ihno3)*xmair/xmhno3/m_dat(region)%data(i,j,1)!mole mole-1
                else if (trim(allvars(j_var)%vname)=='tas'.and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)= temper_dat(1)%data(i,j,1)!K
                else if (trim(allvars(j_var)%vname)=='sfo3'.and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)= tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1)!mole mole-1
!!$                   else if (trim(allvars(j_var)%vname)=='sfpm25'.and. trim(allvars(j_var)%freq) =='1hr')then
!!$                      pm_sizelimit=2.5!micrometers
!!$                      call PMx_Integrate_0d(region,i,j,1,pm_sizelimit,temp,status)!kgm-3
!!$                      allvars(j_var)%data2D(i,j)=  temp/dens/m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfno2'.and. trim(allvars(j_var)%freq) =='1hr')then
                  allvars(j_var)%data2D(i,j)=tracers_c(i,j,1,ino2)*xmair/xmno2/m_dat(region)%data(i,j,1)
                else if (trim(allvars(j_var)%vname)=='sfndtot'.and. trim(allvars(j_var)%freq) =='6hr')then
                  tempbud=0.0
                  dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas)
                  do m=1,7
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens
                    tempbud=tempbud+(tracers(i,j,1,indices(m)))*dens
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens
                  allvars(j_var)%data2D(i,j)=tempbud
                else if (trim(allvars(j_var)%vname)=='sfnd100'.and. trim(allvars(j_var)%freq) =='6hr')then
                  tempbud=0.0
                  dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas)
                  do m=1,7
                    tempbud=tempbud+(1.0-mode_fraction(rw_mode(1,m)%d3(i,j,1),100e-9,m))*(tracers(i,j,1,indices(m)))*dens
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(temp)
                  allvars(j_var)%data2D(i,j)= tempbud
                  !! daily ouptu
                else if (trim(allvars(j_var)%vname)=='sfmmrdust')then
                  tempbud=0.0
                  do m=1,size(allvars(j_var)%indices)
                    index=allvars(j_var)%indices(m)
                    if (index>0)then
                      allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!kg kg-1
                      !tempbud=tempbud+(tracers(i,j,1,index))/m_dat(region)%data(i,j,1)!*den
                    end if
                  end do
                else if (trim(allvars(j_var)%vname)=='depdustabv1')then
                  tempbud=0.0
                  xmcomp=xmdust
                  !modes_dust=(/3,4,6,7/)
                  do m=1,size(allvars(j_var)%indices)
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (allvars(j_var)%indices(m)>0) then
                      if (allvars(j_var)%indices(m)>69)then
                        sedindex=-1
                      else
                        sedindex=sindex(allvars(j_var)%indices(m))
                      end if
                      if (sedindex>0)then
                        do k=1,lmr
                          tempbud=tempbud+( &
                               (deposition(region)%sed(i,j,k,sedindex)*(1.0-mode_fraction(rw_mode(1,modes_dust(m))%d3(i,j,k),1000e-9,modes_dust(m))))) / dxyp(j) * 1.E-03*xmcomp
                        end do
                      end if
                    end if
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud)
                    allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  end do
                else if (trim(allvars(j_var)%vname)=='depdustabv10')then
                  tempbud=0.0
                  xmcomp=xmdust
                  !modes_dust=(/3,4,6,7/)
                  do m=1,size(allvars(j_var)%indices)
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (allvars(j_var)%indices(m)>0) then
                      if (allvars(j_var)%indices(m)>69)then
                        sedindex=-1
                      else
                        sedindex=sindex(allvars(j_var)%indices(m))
                      end if
                      if (sedindex>0)then
                        do k=1,lmr
                          tempbud=tempbud+( &
                               (deposition(region)%sed(i,j,k,sedindex)*(1.0-mode_fraction(rw_mode(1,modes_dust(m))%d3(i,j,k),10000e-9,modes_dust(m))))) / dxyp(j) * 1.E-03*xmcomp
                        end do
                      end if
                    end if
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud)
                    allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  end do
                else if (trim(allvars(j_var)%vname)=='depdustbel1')then
                  tempbud=0.0
                  xmcomp=xmdust
                  modes_dust=(/3,4,6,7/)
                  do m=1,size(allvars(j_var)%indices)
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (allvars(j_var)%indices(m)>0) then
                      if (allvars(j_var)%indices(m)>69)then
                        sedindex=-1
                      else
                        sedindex=sindex(allvars(j_var)%indices(m))
                      end if
                      if (sedindex>0)then
                        do k=1,lmr
                          tempbud=tempbud+( &
                               (deposition(region)%sed(i,j,k,sedindex)*(mode_fraction(rw_mode(1,modes_dust(m))%d3(i,j,k),1000e-9,modes_dust(m))))) / dxyp(j) * 1.E-03*xmcomp
                        end do
                      end if
                    end if
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud)
                    allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  end do
                else if  (trim(allvars(j_var)%vname)=='sfcWind')then
                  yres = dy/yref(1)
                  xres = dx/xref(1)
                  lat = ybeg(1) + 0.5 * yres + yres * (j-1)
                  dxx = ae * xres * gtor * cos(lat*gtor)
                  dyy = ae * yres * gtor
                  vwind= dyy   *(pv_dat(1)%data(i,j,1) + pv_dat(1)%data(i,j+1,1))*0.5 / m_dat(1)%data(i,j,1)
                  uwind=dxx*(pu_dat(1)%data(i,j,1) + pu_dat(1)%data(i-1,j,1))*0.5 / m_dat(1)%data(i,j,1)
                  meanwind=(uwind**2+vwind**2)**(1./2.)
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*meanwind
                else if (trim(allvars(j_var)%vname)=='sfndtot')then
                  tempbud=0.0
                  dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas)
                  do m=1,7
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens
                    tempbud=tempbud+(tracers(i,j,1,indices(m)))*dens
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tracers(i,j,1,indices(m)))*dens
                  allvars(j_var)%data2D(i,j)=tempbud
                else if (trim(allvars(j_var)%vname)=='sfnd100')then
                  tempbud=0.0
                  dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas)
                  do m=1,7
                    tempbud=tempbud+(1.0-mode_fraction(rw_mode(1,m)%d3(i,j,1),100e-9,m))*(tracers(i,j,1,indices(m)))*dens
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(temp)
                  allvars(j_var)%data2D(i,j)= tempbud
                else if (trim(allvars(j_var)%vname)=='sfsh')then
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*humid_dat(1)%data(i,j,1)
                else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_1')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,1)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_2')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,2)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_3')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,3)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfmmraerh2o_4')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*(h2o_mode(region,4)%d3(i,j,1))/m_dat(region)%data(i,j,1)!kg kg-1
                else if (trim(allvars(j_var)%vname)=='tas')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K
                else if (trim(allvars(j_var)%vname)=='sfo3')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1)!mole mole-1
!!$                   else if (trim(allvars(j_var)%vname)=='sfpm25')then
!!$                      pm_sizelimit=2.5!micrometers
!!$                      call PMx_Integrate_0d(region,i,j,1,pm_sizelimit,temp,status)!kg
!!$                      allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j) + dtime*temp/m_dat(region)%data(i,j,1)/dens!kg kg-1
                else if (trim(allvars(j_var)%vname)=='sfno2')then
                  allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j) + dtime*tracers_c(i,j,1,ino2)*xmair/xmno2/m_dat(region)%data(i,j,1)
                else if (trim(allvars(j_var)%vname)=='wetdms' .or.&
                     trim(allvars(j_var)%vname)=='wetno3' .or.&
                     trim(allvars(j_var)%vname)=='wetnh4' .or.&
                     trim(allvars(j_var)%vname)=='wethno3' .or.&
                     trim(allvars(j_var)%vname)=='wetdust' ) then
                  ! Wet index not used for CRESCENDO vars
                  ! Changed to own budget variable to fix budget calculation on the
                  ! variable holder in allvars(jvar). This way you can use the same var
                  ! in more than one time means (hr, daily, monthly)
                  select case ( trim(allvars(j_var)%vname))
                  case('wetdms')
                    !wet_index=6
                    xmcomp=xmdms
                  case('wetno3')
                    !wet_index=11
                    xmcomp=xmno3
                  case('wethno3')
                    !wet_index=12
                    xmcomp=xmhno3
                  case('wetnh4')
                    !wet_index=10
                    xmcomp=xmnh4
                  case('wetdust')
                    !wet_index=5
                    xmcomp=xmdust
                  case DEFAULT
                    write(gol,'("CRESCENDO wetdep xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr
                  end select
                  tempbud=0.0
                  do m=1,n_indices
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (allvars(j_var)%indices(m)<70)then
                      if (allvars(j_var)%indices(m)<=0) then
                        cycle
                      else
                        ! go through gridpoints
                        tempbud=tempbud+(sum(wetdep(region)%lsp(i,j,:,indices(m))) +   sum(wetdep(region)%cp(i,j,:,indices(m)))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                      end if
                    end if
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-wetdepos(region)%f2dslast(i,j,wet_index))
                  !wetdepos(region)%f2dslast(i,j,wet_index)=tempbud
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                  allvars(j_var)%budgetdata(i,j,1)=tempbud
                else if (trim(allvars(j_var)%vname)=='seddust')then
                  tempbud=0.0
                  xmcomp=xmdust
                  do m=1,n_indices
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (indices(m)<=0) then
                      cycle
                    else
                      if (allvars(j_var)%indices(m)>69)then
                        sedindex=-1
                      else
                        sedindex=sindex(allvars(j_var)%indices(m))
                      end if
                      if (sedindex>0)then
                        tempbud=tempbud+( &
                             sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                      end if
                    end if
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud)
                    allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  end do
                else if (trim(allvars(j_var)%vname)=='drydms' .or.&
                     trim(allvars(j_var)%vname)=='dryno3' .or.&
                     trim(allvars(j_var)%vname)=='drynh4' .or.&
                     trim(allvars(j_var)%vname)=='dryhno3' .or.&
                     trim(allvars(j_var)%vname)=='dryno2' .or.&
                     trim(allvars(j_var)%vname)=='dryno' .or.&
                     trim(allvars(j_var)%vname)=='drypan' .or.&
                     trim(allvars(j_var)%vname)=='drydust' &
                     ) then
                  ! Dry index not used for CRESCENDO vars
                  ! Changed to own budget variable to fix budget calculation on the
                  ! variable holder in allvars(jvar). This way you can use the same var
                  ! in more than one time means (hr, daily, monthly)
                  select case ( trim(allvars(j_var)%vname))
                  case('drynh4')
                    !dry_index=10
                    xmcomp=xmnh4
                  case('dryno3')
                    !dry_index=13
                    xmcomp=xmno3
                  case('drydms')
                    !dry_index=6
                    xmcomp=xmdms
                  case('dryhno3')
                    !dry_index=14
                    xmcomp=xmhno3
                  case('dryno2')
                    !dry_index=15
                    xmcomp=xmno2
                  case('dryno')
                    !dry_index=16
                    xmcomp=xmno
                  case('drypan')
                    !dry_index=17
                    xmcomp=xmpan
                  case('drydust')
                    !dry_index=5
                    xmcomp=xmdust
                  case DEFAULT
                    write(gol,'("CRESCENDO drydep xm-var missing: ",a)') trim(allvars(j_var)%vname); call goPr
                  end select
                  tempbud=0.0
                  do m=1,n_indices
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (indices(m)<=0) then
                      cycle
                    else
                      if (allvars(j_var)%indices(m)>69)then
                        sedindex=-1
                      else
                        sedindex=sindex(allvars(j_var)%indices(m))
                      end if
                      !only M7 aerosol tracers (and msa/nh4/no3) are sedimented and deposited
                      if (sedindex>0)then
                        tempbud=tempbud+(buddry_dat(region)%dry(i,j,indices(m)) +&
                             sum(deposition(region)%sed(i,j,:,sedindex))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                      else
                        tempbud=tempbud+buddry_dat(region)%dry(i,j,indices(m)) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                      end if
                    end if
                  end do
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-drydepos(region)%f2dslast(i,j,dry_index))
                  !drydepos(region)%f2dslast(i,j,dry_index)=tempbud
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+(tempbud-allvars(j_var)%budgetdata(i,j,1))
                  allvars(j_var)%budgetdata(i,j,1)=tempbud
                else if (trim(allvars(j_var)%vname)=='emioa')then
                  emi_index=1
                  xmcomp=xmpom
                  tempbud=0.0
                  do m=1,n_indices
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (indices(m)<=0) then
                      cycle
                    else
                      tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                      !OLDWAY which is not in use for CRESCENDO vars:
                      !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-emission(region)%f2dslast(i,j,emi_index))
                      !emission(region)%f2dslast(i,j,emi_index)=tempbud
                    end if
                  end do
                  ! emioa is requested as emi+chepsoa so add soa production
                  !if (trim(allvars(j_var)%freq)=='mon')then
                  !   tempbud= tempbud + AC_diag_prod(region)%prod(i,j,k,3)/ dxyp(j) !kg m-2s-1 (s-1 at the time of writing)
                  !else if (trim(allvars(j_var)%freq)=='6hr')then
                  !   tempbud= tempbud + AC_diag_prod(region)%prod(i,j,k,4)/ dxyp(j) !kg m-2s-1 (s-1 at the time of writing)
                  !else
                  !   write(gol,'("Uknown output frequency for production of soa in calculation of : ",a)')trim(allvars(j_var)%freq); call goErr
                  !   TRACEBACK
                  !   status=1; return
                  !end if
                  if (trim(allvars(j_var)%freq)=='6hr'.or.trim(allvars(j_var)%freq)=='1hr')then
                    allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep   /dtime
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then
                    allvars(j_var)%data2D(i,j)=  allvars(j_var)%data2D(i,j)+tempbud-allvars(j_var)%budgetdata(i,j,1)
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  else
                    write(gol,'("Uknown output frequency: ",a)')trim(allvars(j_var)%vname); call goErr
                    TRACEBACK
                    status=1; return
                  end if
                else if (trim(allvars(j_var)%vname)=='emiisop')then
                  tempbud=0.0
                  gasind=iisop
                  xmcomp=xmisop
                  emi_index=8
                  tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp !kg m-2s-1 (s-1 at the time of writing)
                  if (trim(allvars(j_var)%freq)=='6hr'.or.trim(allvars(j_var)%freq)=='1hr')then
                    allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime! instantaneous values -> /dtime
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing)
                  else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then
                    allvars(j_var)%data2D(i,j)=  allvars(j_var)%data2D(i,j) + (tempbud-allvars(j_var)%budgetdata(i,j,1))
                    !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing)
                  else
                    write(gol,'("Uknown output frequency: ",a)')trim(allvars(j_var)%vname); call goErr
                    TRACEBACK
                    status=1; return
                  end if
                  !allvars(j_var)%data2D(i,j)= tempbud-allvars(j_var)%budgetdata(i,j,1)
                  !allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(tempbud-allvars(j_var)%budgetdata(i,j,1))
                  !allvars(j_var)%budgetdata(i,j,1)=tempbud
                else if (trim(allvars(j_var)%vname)=='emimono')then
                  tempbud=0.0
                  gasind=iterp
                  xmcomp=xmterp
                  emi_index=13
                  tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp
                  if (trim(allvars(j_var)%freq)=='6hr'.or.trim(allvars(j_var)%freq)=='1hr')then
                    !Instantaneous values requested
                    allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime! instantaneous values -> /dtime : kgm-2s-1
                    allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing)
                  else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then
                    !mean values
                    allvars(j_var)%data2D(i,j)=  allvars(j_var)%data2D(i,j) + (tempbud-allvars(j_var)%budgetdata(i,j,1))
                    allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing)
                  else
                    write(gol,'("Uknown output frequency: ",a)')trim(allvars(j_var)%vname); call goErr
                    TRACEBACK
                    status=1; return
                  end if
                else if (trim(allvars(j_var)%vname)=='emidms')then
                  tempbud=0.0
                  gasind=idms
                  emi_index=6
                  xmcomp=xmdms
                  tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,gasind)) / dxyp(j) * 1.E-03*xmcomp
                  if (trim(allvars(j_var)%freq)=='6hr')then
                    allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime
                    allvars(j_var)%budgetdata(i,j,1)=tempbud !kg m-2s-1 (s-1 at the time of writing)
                  else
                    write(gol,'("CRESCENDO output frequency not implemented: ",a)')trim(allvars(j_var)%vname); call goErr
                    TRACEBACK
                    status=1; return
                  end if
                else if (trim(allvars(j_var)%vname)=='emiss')then
                  emi_index=7
                  xmcomp=xmnacl
                  do m=1,n_indices
                    !exit if indices ==0
                    !since indices after 0 will all be also 0
                    if (indices(m)<=0) then
                      cycle
                    else
                      ! go through gridpoints
                      tempbud=tempbud+ sum(budemi_dat(region)%emi(i,j,:,indices(m))) / dxyp(j) * 1.E-03*xmcomp
                    end if
                  end do
                  if (trim(allvars(j_var)%freq)=='6hr')then
                    allvars(j_var)%data2D(i,j)= (tempbud-allvars(j_var)%budgetdata(i,j,1))/dtime!timerange for instantaneous data is half of dynamical timestep /dtime !kg m-2s-1 (s-1 at the time of writing)
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  else if (trim(allvars(j_var)%freq)=='mon'.or.trim(allvars(j_var)%freq)=='day')then
                    allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j) + (tempbud-allvars(j_var)%budgetdata(i,j,1))*dtime
                    allvars(j_var)%budgetdata(i,j,1)=tempbud
                  end if
                else if (trim(allvars(j_var)%unit)=='cm-3') then
                  dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas)
                  !1e-6 for m3->cm-3
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,index)*dens*1e-6
                else if (index>0 .and. index<=ntracet) then
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,index)/m_dat(region)%data(i,j,1)!kg kg-1
                else if (index>ntracet) then ! than non transported (tracers_c=> chem_dat(region)%rm)
                  allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*tracers_c(i,j,1,index)/m_dat(region)%data(i,j,1)!kg kg-1
                else
                  write(gol,*) 'user_output_aerchemmip: missing accumulation for CRESCENDO 3D variable ',allvars(j_var)%vname ; call goErr
                  TRACEBACK
                  status=1; return
                end if
              end do
            end do
          end if

          !end if
        else if (trim(allvars(j_var)%class)=='mmr')then
          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)
                !do all mode-specific tracers for given compound
                dens = pres3d(i,j,k) / temper_dat(region)%data(i,j,k) * xmair * 1.E-3 / (m_dat(region)%data(i,j,k) * rgas)
                do imode=1,size(allvars(j_var)%indices)
                  index=allvars(j_var)%indices(imode)
                  ! first variables that that are not tracers (transported or non-transported)
                  if (index<=0) then
                    if (trim(allvars(j_var)%vname)=='mmrpm2p5')then
                      pm_sizelimit=2.5!micrometers
                      !call PMx_Integrate_3d(region,pm_sizelimit,temp_pm,status)
                      CALL pmx_integrate_0d( region, i, j, k, pm_sizelimit, temp, status ) !returns kg/m3
                      allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temp/dens/m_dat(region)%data(i,j,k) !kg(pm)/m3 /(1/m3) /kg(air)->kg(pm)/kg(air)
                    else if (trim(allvars(j_var)%vname)=='mmrpm1')then
                      pm_sizelimit=1.0!micrometers
                      !call PMx_Integrate_3d(region,pm_sizelimit,temp_pm,status)
                      CALL pmx_integrate_0d( region, i, j, k, pm_sizelimit, temp, status ) !returns kg/m3
                      allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temp/dens/m_dat(region)%data(i,j,k)!kg kg-1
                    else if (trim(allvars(j_var)%vname)=='mmrpm10')then
                      pm_sizelimit=10.0!micrometers
                      !call PMx_Integrate_3d(region,pm_sizelimit,temp_pm,status)
                      CALL pmx_integrate_0d( region, i, j, k, pm_sizelimit, temp, status ) !returns kg/m3
                      !kg/m3/m-3/kg
                      allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temp/dens/m_dat(region)%data(i,j,k)!kg kg-1
                    else if (trim(allvars(j_var)%vname)=='mmraerh2o')then
                      allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*(h2o_mode(region,1)%d3(i,j,k)+h2o_mode(region,2)%d3(i,j,k)+h2o_mode(region,3)%d3(i,j,k)+h2o_mode(region,4)%d3(i,j,k))/m_dat(region)%data(i,j,k)!kg kg-1
                    end if
                    !exit
                  else
                    !Transported variables are stored in different array (tracers=> mass_dat(region)%rm)
                    if (index<=ntracet) then
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1
                    else ! than non transported (tracers_c=> chem_dat(region)%rm)
                      allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)/m_dat(region)%data(i,j,k)!kg kg-1
                    end if
                  end if
                end do
              end do
            end do
          end do

          !end if
        else if (trim(allvars(j_var)%class)=='gas')then
          do k=1,lmr
            do j=j1,j2
              do i=i1,i2
                index=allvars(j_var)%indices(1)
                if (index<=0) then
                  cycle
                else  if (trim(allvars(j_var)%vname)=='h2o') then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*humid_dat(1)%data(i,j,k)/xmh2o*xmair !mole mole-1
                  !tracers_c(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k)
                else if (trim(allvars(j_var)%vname)=='voc')then
                  ! ACTUALLY not requested
                  do m=1,ncb5
                    gasind=nmhc(m)
                    xmgas=xmcb5(m)
                    select case(nmhc(m))
                    case(ipar,ich2o,ich3oh,ihcooh)
                      xmcomp=xmc
                    case(ieth,iole,iald2,imcooh,ic2h6,iethoh)
                      xmcomp=xmc*2.0
                    case(imgly,ic3h8,ic3h6,iacet)
                      xmcomp=xmc*3.0
                    end select
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,gasind)*xmair/xmcomp/m_dat(region)%data(i,j,k) !kmole kmole-1
                  end do
                  do m=1,2
                    if (m==1)then
                      gasind=iterp
                      xmcomp=xmc*5*2
                    else
                      gasind=iisop
                      xmcomp=xmc*5.0
                    end if
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,gasind)*xmair/xmcomp/m_dat(region)%data(i,j,k) !kmole kmole-1
                  end do
                else
                  xmgas=allvars(j_var)%xmgas
                  if (index >= ntracet) then
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k) !kmole kmole-1
                  else
                    allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)*xmair/xmgas/m_dat(region)%data(i,j,k) !kmole kmole-1
                  end if
                end if
              end do
            end do
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='ps6h')then
          do i=i1,i2
            do j=j1,j2
              if (trim(allvars(j_var)%vname)=='ps')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa
              end if
            end do
          end do
          !end if

          ! 1 hourly surface variables
        else if (trim(allvars(j_var)%class)=='sf1h')then
          do i=i1,i2
            do j=j1,j2
              if (trim(allvars(j_var)%vname)=='ps')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa
              else if (trim(allvars(j_var)%vname)=='sfno2')then
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*tracers_c(i,j,1,ino2)*xmair/xmno2/m_dat(region)%data(i,j,1)*dtime !kmole kmole-1
              else if (trim(allvars(j_var)%vname)=='sfo3')then
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1) !kmole kmole-1
              else if (trim(allvars(j_var)%vname)=='sfpm25')then
                pm_sizelimit=2.5!micrometers
                dens = pres3d(i,j,1) / temper_dat(region)%data(i,j,1) * xmair * 1.E-3 / (m_dat(region)%data(i,j,1) * rgas)
                call PMx_Integrate_0d(region,i,j,1,pm_sizelimit,temp,status)
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temp/m_dat(region)%data(i,j,1)/dens!kg kg-1
                !allvars(j_var)%data2D(i,j)=temp/m_dat(region)%data(i,j,1)/dens*dtime!kg kg-1
              else if (trim(allvars(j_var)%vname)=='tas')then
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K
              end if
            end do
          end do

          !end if
          ! surface daily variables
        else if (trim(allvars(j_var)%class)=='sf1d')then
          do i=i1,i2
            do j=j1,j2
              if (trim(allvars(j_var)%vname)=='ps')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*sp_dat(1)%data(i,j,1)!Pa
              else if (trim(allvars(j_var)%vname)=='toz')then
                ! conversion (in order of execution):
                !kg->g (1e3), g->molec (xmo3), m2->cm2(1e-4) , molec-> molec/m2(dxyp) , molec/cm2->dobson (DOBS)
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*sum(tracers(i,j,:,io3)*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs) !GU
              else if (trim(allvars(j_var)%vname)=='sfo3max')then
                sfo3_molmol = tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1)
                if ( sfo3_molmol > allvars(j_var)%data2D(i,j)) then
                  allvars(j_var)%data2D(i,j) = sfo3_molmol !kmole kmole-1
                end if
              else if (trim(allvars(j_var)%vname)=='maxpblz')then
                if (conv_dat(1)%blh(i,j)>  allvars(j_var)%data2D(i,j)) then
                  allvars(j_var)%data2D(i,j)= conv_dat(1)%blh(i,j)!m
                end if
              else if (trim(allvars(j_var)%vname)=='minpblz')then
                if (conv_dat(1)%blh(i,j)<  allvars(j_var)%data2D(i,j)) then
                  allvars(j_var)%data2D(i,j)= conv_dat(1)%blh(i,j)!m
                end if
              else if (trim(allvars(j_var)%vname)=='tas')then
                allvars(j_var)%data2D(i,j)=allvars(j_var)%data2D(i,j)+dtime*temper_dat(1)%data(i,j,1)!K
              else if (trim(allvars(j_var)%vname)=='tasmin')then
                if (temper_dat(1)%data(i,j,1)<  allvars(j_var)%data2D(i,j)) then
                  allvars(j_var)%data2D(i,j)= temper_dat(1)%data(i,j,1)!K
                end if
              else if (trim(allvars(j_var)%vname)=='tasmax')then
                if (temper_dat(1)%data(i,j,1)>  allvars(j_var)%data2D(i,j)) then
                  allvars(j_var)%data2D(i,j)= temper_dat(1)%data(i,j,1)!K
                end if
              else if  (trim(allvars(j_var)%vname)=='pr')then
                allvars(j_var)%data2D(i,j)= allvars(j_var)%data2D(i,j)+dtime*(lsp_dat(1)%data(i,j,1)+ cp_dat(1)%data(i,j,1))/rol!kgm-2s-1
              end if
            end do
          end do
          !end if

        else if (trim(allvars(j_var)%class)=='zonal')then
          ! zonal mean needs to be done afterwards...
          ! So just save the variables as 3D
          do i=i1,i2
            do j=j1,j2
              do k=1,lmr
                if (trim(allvars(j_var)%vname)=='ch4')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,ich4)/m_dat(region)%data(i,j,k)*xmair/xmch4!mol mol-1
                else if (trim(allvars(j_var)%vname)=='hno3')then
                  allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,ihno3)/m_dat(region)%data(i,j,k)*xmair/xmhno3!mol mol-1
                else if (trim(allvars(j_var)%vname)=='o3')then
                  allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,io3)/m_dat(region)%data(i,j,k)*xmair/xmo3!mol mol-1
                else if (trim(allvars(j_var)%vname)=='ta')then
                  allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*temper_dat(1)%data(i,j,k)!K
                else if (trim(allvars(j_var)%vname)=='h2o')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+0.0!dtime*humid_dat(1)%data(i,j,k)/xmh2o*xmair!mol mol-1
                else if (trim(allvars(j_var)%vname)=='zg')then
                  allvars(j_var)%data3D(i,j,k)= allvars(j_var)%data3D(i,j,k)+dtime*gph_dat(1)%data(i,j,k)!m
                else if (trim(allvars(j_var)%vname)=='ho2')then
                  allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,iho2)/m_dat(region)%data(i,j,k)*xmair/xmho2!mol mol-1
                else if (trim(allvars(j_var)%vname)=='oh')then
                  allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,ioh)/m_dat(region)%data(i,j,k)*xmair/xmoh!mol mol-1
                else if (trim(allvars(j_var)%vname)=='noy')then
                  do icomp=1,size( allvars(j_var)%indices(:))
                    index= allvars(j_var)%indices(icomp)
                    xmcomp=ra(index)
                    if (allvars(j_var)%indices(icomp)>0) then
                      !write(gol,'(": ",a,i3,E20.7e2)') trim(allvars(j_var)%vname),index,tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp; call goPr
                      if  (allvars(j_var)%indices(icomp)<70) then
                        allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp!???
                      else
                        allvars(j_var)%data3D(i,j,k)=allvars(j_var)%data3D(i,j,k)+dtime*tracers_c(i,j,k,index)/m_dat(region)%data(i,j,k)*xmair/xmcomp!???
                      end if
                    end if
                  end do
                end if
              end do
            end do
          end do

        else
          ! optics and fixed are not accumulated here
          ! optics in optics_calc
          ! fixed not at all
          if (trim(allvars(j_var)%class)/='optics' .and. trim(allvars(j_var)%class)/='fixed') then
            write(gol,*) 'output_aerhemmip_accumulate: output class not found!!!',trim(allvars(j_var)%vname),trim(allvars(j_var)%class)
            !call goPr
            call goErr
            status=1
            return
          end if
        end if
      end do
    end if

    ! zero accumulated budget variables for the amount between output steps
    AC_diag_prod(region)%prod(i1:i2,j1:j2,:,1)=0.0
    AC_diag_prod(region)%prod(i1:i2,j1:j2,:,2)=0.0
    AC_diag_prod(region)%prod(i1:i2,j1:j2,:,3)=0.0
    AC_diag_prod(region)%prod(i1:i2,j1:j2,:,4)=0.0
    AC_loss(region)%prod(i1:i2,j1:j2,:,1:2)=0.0
    AC_O3_lp(region)%prod(i1:i2,j1:j2,:,1:2)=0.0
    deallocate(gphx)
    deallocate(tx)
    deallocate(temp_pm)
    deallocate(pres3d)
    deallocate(pres3dh)

    call GO_Timer_Start( itim_accu_opt, status )

    call calculate_optics(1,dhour,l_do_ec550aer_only,status)
    call GO_Timer_End( itim_accu_opt, status )
    call GO_Timer_End( itim_accu, status )
    !    IF_NOTOK_RETURN(status=1)
    call goLabel()
    !status = 1
  end subroutine accumulate_data

  subroutine output_aerchemmip_done(status)
    use partools, only: isRoot,myid
    implicit none
    integer :: status
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_done'
    integer :: j_var, region

    call goLabel(rname)
    region = 1
    if (isroot) then
      do j_var=1,n_vars
        call MDF_Close( allvars(j_var)%fileunit, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      end do
    end if

    deallocate(wdep_out )
    deallocate(drydepos(region)%f2dslast) 
    deallocate(wetdepos(region)%f2dslast) 
    deallocate(emission(region)%f2dslast) 
    deallocate(dimension_data%lon) 
    deallocate(dimension_data%lat) 
    deallocate(dimension_data%lev) 

    do j_var=1,n_vars
      deallocate(allvars(j_var)%data2D)
      deallocate(allvars(j_var)%data3D)
    end do
    deallocate(allvars)
    deallocate(fixedvars)

    call goLabel()
    status = 0
  end subroutine output_aerchemmip_done

  !subroutine add_variable(itm5,shortname,longname,unit,positive,standardname,varid,fileunit,filename,data_dims,status,class,table,pxmgas)
  subroutine add_variable(itm5,shortname,longname,unit,data_dims,status,class,table,pxmgas)
    use chem_param,  only: mode_end_so4,mode_end_pom,mode_end_bc,mode_end_ss,mode_end_dust
    use global_data, only: outdir
    use datetime,    only: tau2date, date2tau
    use dims,        only: itau,itaue,itaut

    implicit none
    integer:: itm5
    character(len=*),intent(in)::shortname
    character(len=*),intent(in)::longname
    character(len=*)::unit
    character(len=30)::standardname

    character(len=*)::table
    character(len=*),optional::class
    real,optional::pxmgas
    integer:: data_dims
    integer,intent(out)::status
    !LOCAL
    character(len=4)::positive=''
    integer:: fileunit=-1 !defined only when creating a file
    integer:: varid=-1! defined when opening ncfile
    !character(len=120)::filename
    character(len=30)::table_id
    !character(len=30)::source_id
    !character(len=30)::experiment_id
    character(len=30)::member_id
    !character(len=30)::grid_label
    character(len=30)::time_range
    character(len=200)::filename1
    character(len=10)::freq
    real,dimension(:,:),pointer::data2D
    !    real,dimension(:,:),pointer::dataZonal
    real,dimension(:,:,:),pointer::data3D
    real,dimension(:,:,:),pointer::budget

    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_add_variable'
    integer ::i1,i2,j1,j2,jmr,imr,lmr
    integer, dimension(6) :: idater, idateend, idatett
    integer               :: endmonth, endday
    character(len=30)     :: idates

    call tau2date(itau,idater)

    ! define frequency from table
    if (trim(table)=='AERhr')then
      !table id depends on variable
      table_id=table
      freq='1hr'
    else if (trim(table)=='AER6hr')then
      !table id depends on variable
      table_id=table
      freq='6hr'
    else if( trim(table)=='AERmon'.or.trim(table)=='AERmonZ'.or.trim(table)=='Emon')then
      !table id depends on variable
      table_id=table
      freq='mon'
    else if(trim(table)=='AERday')then
      !table id depends on variable
      table_id=table
      freq='day'
    else if(trim(table)=='AERfx')then
      !table id depends on variable
      table_id=table
      freq='na'
    else
      freq='freq-nondefined'
      table_id='table-nondefined'
    end if

    ! CREATE date string for output
    ! 
    ! ATM assumed that the output is initilised at the beginninh of year
    endmonth=12
    endday=31
    !
    if (freq=='mon')then
      ! By default AERCHEMMIP runs are done by 1-year chunks -> idater(2) -  idater(2)+11
      write(idates,  '(i4,i2.2,a,i4,i2.2)') idater(1), idater(2),'-', idater(1),endmonth
    else if(freq=='day')then
      !time range form Jan-1 ->Dec-31x
      write(idates,  '(i4,2i2.2,a,i4,2i2.2)') idater(1), idater(2), idater(3),'-', idater(1), endmonth, endday

    else if(freq=='1hr')then
      write(idates,  '(i4,2i2.2,2a2,a,i4,2i2.2,2a2)') idater(1), idater(2), idater(3),'00','00','-', idater(1), endmonth, endday, '23', '59'    
    else if (freq=='6hr')then

      write(idates,  '(i4,2i2.2,2a2,a,i4,2i2.2,2a2)') idater(1), idater(2), idater(3),'00','00','-', idater(1), endmonth, endday,'18','00'    
    end if
    !statf(region)%fname = trim(outdir)//'/'//&
    !           trim(f_dataid) //'_'//&
    !           trim(f_model)  //'_'//&
    !           trim(aerocom_exper) //'_'//&
    !           trim(f_dimstat)//'_'//&
    !           trim(idates)   //'_'//&
    !           trim(aerocom_freq)   //'.nc'

    call goLabel(rname)

    call GO_Timer_Start( itim_addvar, status )
    IF_NOTOK_RETURN(status=1)
    !outdir='output'


    ! temporary
    standardname=longname
    ! source_id constant
    !source_id='EC-EARTH-AerChem'
    ! experiment depends on run
    !experiment_id='exp_dynamic'
    member_id='r'//trim(realization_i)//'i'//trim(initialization_i)//'p'//trim(physics_i)//'f'//trim(forcing_i)
    !grid_label='3x2_degrees'

    ! time range has divider in place since it can be omitted alltogether with non-time dependendent variables
    if (trim(table)=='AERfx')then
      time_range=''
    else
      time_range='_'//trim(idates)
    end if

    ! for fixed variables time range should not be written
    n_vars=n_vars+1
    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
    if (trim(shortname)=='phalf') then
      allocate(budget(i1:i2,j1:j2,1:lmr+1))
    else
      allocate(budget(i1:i2,j1:j2,1:lmr))
    end if
    budget(:,:,:)=0.0
    ! 2D variables
    if (data_dims==2) then
      !Allocate holders for data
      allocate(allvars(n_vars)%data2D(i1:i2,j1:j2))
      allocate(allvars(n_vars)%data3D(1,1,1))

      ! allocate local variables
      allocate(data2D(i1:i2,j1:j2))
      allocate(data3D(1,1,1))
      ! zero local data holders
      data2D(:,:)=0.0
      !       dataZonal(:,:)=0.0
      data3D(:,:,:)=0.0
      !init variable

      ! set default for minimum variables to high value
      if (shortname=='minpblz' .or. shortname=='tasmin')then
        data2D(:,:)=1000000.0
      end if
      !create filename

      if (trim(class)=='crescendo')then

        !          filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//'_'//trim(time_range)//trim('.nc')
        filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(class)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')
      else
        filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')
      end if
      allvars(n_vars)=varfile(itm5,shortname,longname,unit,positive,standardname,data2D,data3D,budget,-1,-1,&
           filename1,2,-1,-1,-1,-1,-1,-1,-1,-1,-1,data_dims,freq,class,(/0,0,0,0,0,0,0,0,0,0,0/),pxmgas,table_id)


      !! LEFT HERE on purpose to see what variables go where in above statement
!!$       allvars(n_vars)%itm5=itm5
!!$       allvars(n_vars)%vname=shortname
!!$       allvars(n_vars)%lname=longname
!!$       allvars(n_vars)%unit=unit
!!$       allvars(n_vars)%positive=positive
!!$       allvars(n_vars)%standard_name=standardname
!!$       allvars(n_vars)%data2D=data2D
!!$       allvars(n_vars)%data3D=data3D
!!$       allvars(n_vars)%budgetdata=data3D
!!$       allvars(n_vars)varid=-1
!!$       allvars(n_vars)%filenunit=-1
!!$       allvars(n_vars)%filename=filename1
!!$       allvars(n_vars)%dimensions=2
!!$       allvars(n_vars)%lon_varid=-1
!!$       allvars(n_vars)%lat_varid=-1
!!$       allvars(n_vars)%lev_varid=-1
!!$       allvars(n_vars)%time_varid=-1
!!$       allvars(n_vars)%bounds_varid=-1
!!$       allvars(n_vars)%dims=dims
!!$       allvars(n_vars)%class=class
!!$       allvars(n_vars)%indices=(/0,0,0,0,0,0,0/)) 
!!$       allvars(n_vars)%xmgas=molarweight
!!$       allvars(n_vars)%table_id=


      ! 3D variables
    else if (data_dims==3) then
      ! allocate holders for data
      allocate(allvars(n_vars)%data2D(1,1))

      if (trim(shortname)=='phalf') then
        allocate(allvars(n_vars)%data3D(i1:i2,j1:j2,1:lmr+1))
        allocate(data3D(i1:i2,j1:j2,1:lmr+1))

      else
        allocate(allvars(n_vars)%data3D(i1:i2,j1:j2,1:lmr))
        allocate(data3D(i1:i2,j1:j2,1:lmr))


      end if


      ! allocate local variables
      ! maybe remove these
      allocate(data2D(1,1))
      !allocate(data3D(i1:i2,j1:j2,1:lmr))
      ! zero local data holders
      data2D(:,:)=0.0
      data3D(:,:,:)=0.0

      !init variable
      if (trim(class)=='crescendo')then

        !          filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//'_'//trim(time_range)//trim('.nc')
        filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(class)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')
      else
        filename1= trim(outdir)//'/'//trim(shortname)//'_'//trim(table_id)//'_'//trim(source_id)//'_'//trim(experiment_id)//'_'//trim(member_id)//'_'//trim(grid_label)//trim(time_range)//trim('.nc')


      end if
      allvars(n_vars)=varfile(itm5,shortname,longname,unit,positive,standardname,data2D,data3D,budget,-1,-1,&
           filename1,3,-1,-1,-1,-1,-1,-1,-1,-1,-1,data_dims,freq,class,(/0,0,0,0,0,0,0,0,0,0,0/),pxmgas,table)
    end if
    ! add chemical info also: (vars beginning with emi,wet,dry)
    select case (trim(shortname(4:)))
    case ('so4')
      allvars(n_vars)%indices(1:7)=(/iso4nus,iso4ais,iso4acs,iso4cos,0,0,0/)!mode_end_so4
    case ('so2')
      allvars(n_vars)%indices(1)=iso2
    case ('oa','emioa')
      !allvars(n_vars)%indices(1:7)=(/ipomais,ipomacs,ipomcos,ipomaii,0,0,0/)!mode_end_pom
      allvars(n_vars)%indices(1:9)=(/ipomais,ipomacs,ipomcos,ipomaii,isoanus,isoaais,isoaacs,isoacos,isoaaii/)!mode_end_pom
    case ('poa','emipoa')
      allvars(n_vars)%indices(1:4)=(/ipomais,ipomacs,ipomcos,ipomaii/)!mode_end_pom
    case ('soa')
      allvars(n_vars)%indices(1:7)=(/isoanus,isoaais,isoaacs,isoacos,isoaaii,0,0/)!mode_end_pom
    case ('bc')
      allvars(n_vars)%indices(1:7)=(/ibcaii,ibcais,ibcacs,ibccos,0,0,0/)!mode_end_bc

    case ('ss','emiss')
      allvars(n_vars)%indices(1:7)=(/issacs,isscos,0,0,0,0,0/)!mode_end_ss
    case ('dust')
      allvars(n_vars)%indices(1:7)=(/iduacs,iducos,iduaci,iducoi,0,0,0/)!mode_end_dust
    case('nox')
      !allvars(n_vars)%indices(1:2)=(/ino,ino2/)
      allvars(n_vars)%indices(1)=inox
    case('voc')
      allvars(n_vars)%indices(1)=1!(/ipar,ieth,iole,iald2,imgly,0,0/)
    case('bvoc')
      allvars(n_vars)%indices(1:2)=(/iterp,iisop/)!ibvoc
    case('nh3','sfnh3')
      allvars(n_vars)%indices(1)=inh3
    case('nh4')
      allvars(n_vars)%indices(1)=inh4
    case('noy')
      allvars(n_vars)%indices(1:11)=(/ ino, iNO2, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/)
      !allvars(n_vars)%indices(1)=ino2!(/ ino, iNO2, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/)
    case('no3')
      allvars(n_vars)%indices(1)=ino3_a
    case('pm1')
      allvars(n_vars)%indices(1)=-1
    case('pm2p5','sfpm25')
      allvars(n_vars)%indices(1)=-1
    case('pm10')
      allvars(n_vars)%indices(1)=-1
    case('o3')
      allvars(n_vars)%indices(1)=io3
    case('co')
      allvars(n_vars)%indices(1)=ico
    case('dms')
      allvars(n_vars)%indices(1)=idms
    case('isop')
      allvars(n_vars)%indices(1)=iisop
      !case('terp')
      !   allvars(n_vars)%indices(1)=iterp

    end select
    select case (trim(shortname))
    case('noy')
      allvars(n_vars)%indices(1:11)=(/ iNO, ino2, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/)
      !allvars(n_vars)%indices(1)=ino2!(/ iNO, ino3, iHNO3, iNO3_a, ihno4, in2o5, iPAN, iOrgNtr, ihono, ich3o2no2/)
    case('areacella')
      allvars(n_vars)%indices(:)=0
      areacella=n_vars

    case('c2h2')
      allvars(n_vars)%indices(1)=-1
    case('c2h6')
      allvars(n_vars)%indices(1)=ic2h6
    case('c3h6')
      allvars(n_vars)%indices(1)=ic3h6
    case('c3h8')
      allvars(n_vars)%indices(1)=ic3h8
    case('ch3coch3')
      allvars(n_vars)%indices(1)=iacet!ich3coch3
    case('ch4')
      allvars(n_vars)%indices(1)=ich4
    case('co')
      allvars(n_vars)%indices(1)=ico
    case('so2')
      allvars(n_vars)%indices(1)=iso2
    case('co2')
      allvars(n_vars)%indices(1)=-1
    case('dms')
      allvars(n_vars)%indices(1)=idms
    case('h2o')
      allvars(n_vars)%indices(1)=-1!ih2o
    case('hcho')
      allvars(n_vars)%indices(1)=ich2o
    case('hcl')
      allvars(n_vars)%indices(1)=-1!ihcl
    case('hno3','sfmmrhno3')
      allvars(n_vars)%indices(1)=ihno3
    case('isop')
      allvars(n_vars)%indices(1)=iisop
    case('n2o')
      allvars(n_vars)%indices(1)=-1!in2o
    case('no', 'sfno')
      allvars(n_vars)%indices(1)=ino
    case('no2','sfno2')
      allvars(n_vars)%indices(1)=ino2
    case('o3','sfo3')
      allvars(n_vars)%indices(1)=io3
    case('o3ste')
      allvars(n_vars)%indices(1)=io3s
    case('oh')
      allvars(n_vars)%indices(1)=ioh
    case('pan')
      allvars(n_vars)%indices(1)=ipan
      !CRESCENDO variables
    case('nh3','sfmmrnh3','sfnh3')
      allvars(n_vars)%indices(1)=inh3
    case('nh4','sfmmrnh4')
      allvars(n_vars)%indices(1)=inh4
    case('sfmmrno3')
      allvars(n_vars)%indices(1)=ino3_a
    case('pm1')
      allvars(n_vars)%indices(1)=-1
    case('pm2p5','sfpm25')
      allvars(n_vars)%indices(1)=-1
    case ('mmraerh2o_1','mmraerh2o_2','mmraerh2o_3','mmraerh2o_4','mmraerh2o')
      allvars(n_vars)%indices(1)=-1
    case ('sfmmrso4')
      allvars(n_vars)%indices(1:4)=(/iso4nus,iso4ais,iso4acs,iso4cos/)
    case ('sfmmrsoa')
      allvars(n_vars)%indices(1:5)=(/isoanus,isoaais,isoaacs,isoacos,isoaaii/)
    case ('sfmmroa')
      allvars(n_vars)%indices(1:4)=(/ipomais,ipomacs,ipomcos,ipomaii/)
    case ('sfmmrbc')
      allvars(n_vars)%indices(1:4)=(/ibcais,ibcacs,ibccos,ibcaii/)
    case ('sfmmrdustbel1','sfmmrdustabv1','sfmmrdustabv10','seddust','sfmmrdust')
      allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/)
    case ('sfmmrss')
      allvars(n_vars)%indices(1:2)=(/issacs,isscos/)
    case ('mmrsu1','sfmmrsu1')
      allvars(n_vars)%indices(1)=iso4nus
    case ('mmrsu2','sfmmrsu2')
      allvars(n_vars)%indices(1)=iso4ais
    case ('mmrsu3','sfmmrsu3')
      allvars(n_vars)%indices(1)=iso4acs
    case ('mmrsu4','sfmmrsu4')
      allvars(n_vars)%indices(1)=iso4cos
    case ('mmrsoa1','sfmmrsoa1')
      allvars(n_vars)%indices(1)=isoanus
    case ('mmrsoa2','sfmmrsoa2')
      allvars(n_vars)%indices(1)=isoaais
    case ('mmrsoa3','sfmmrsoa3')
      allvars(n_vars)%indices(1)=isoaacs
    case ('mmrsoa4','sfmmrsoa4')
      allvars(n_vars)%indices(1)=isoacos
    case ('mmrsoa5','sfmmrsoa5')
      allvars(n_vars)%indices(1)=isoaaii
    case ('mmroa2','sfmmroa2')
      allvars(n_vars)%indices(1)=ipomais
    case ('mmroa3','sfmmroa3')
      allvars(n_vars)%indices(1)=ipomacs
    case ('mmroa4','sfmmroa4')
      allvars(n_vars)%indices(1)=ipomcos
    case ('mmroa5','sfmmroa5')
      allvars(n_vars)%indices(1)=ipomaii
    case ('mmrbc2','sfmmrbc2')
      allvars(n_vars)%indices(1)=ibcais
    case ('mmrbc3','sfmmrbc3')
      allvars(n_vars)%indices(1)=ibcacs
    case ('mmrbc4','sfmmrbc4')
      allvars(n_vars)%indices(1)=ibccos
    case ('mmrbc5','sfmmrbc5')
      allvars(n_vars)%indices(1)=ibcaii
    case ('mmrss3','sfmmrss3')
      allvars(n_vars)%indices(1)=issacs
    case ('mmrss4','sfmmrss4')
      allvars(n_vars)%indices(1)=isscos
    case ('mmrdu3','sfmmrdu3')
      allvars(n_vars)%indices(1)=iduacs
    case ('mmrdu4','sfmmrdu4')
      allvars(n_vars)%indices(1)=iducos
    case ('mmrdu6','sfmmrdu6')
      allvars(n_vars)%indices(1)=iduaci
    case ('mmrdu7','sfmmrdu7')
      allvars(n_vars)%indices(1)=iducoi
    case ('sfndtot','sfnd100','ndtot')
      allvars(n_vars)%indices(1:7)=(/inus_n,iais_n,iacs_n,icos_n,iaii_n,iaci_n,icoi_n/)
    case ('nd1','sfnd1','sfnd1_tst')
      allvars(n_vars)%indices(1)=inus_n
    case ('nd2','sfnd2')
      allvars(n_vars)%indices(1)=iais_n
    case ('nd3','sfnd3')
      allvars(n_vars)%indices(1)=iacs_n
    case ('nd4','sfnd4')
      allvars(n_vars)%indices(1)=icos_n
    case ('nd5','sfnd5')
      allvars(n_vars)%indices(1)=iaii_n
    case ('nd6','sfnd6')
      allvars(n_vars)%indices(1)=iaci_n
    case ('nd7','sfnd7')
      allvars(n_vars)%indices(1)=icoi_n
      !  case('sfmmrnh3')
      !    allvars(n_vars)%indices(1)=inh3
      !    case('sfmmrnh4')
      !       allvars(n_vars)%indices(1)=inh4
    case('drydms')
      allvars(n_vars)%indices(1)=idms
    case('wetdms')
      allvars(n_vars)%indices(1)=idms
    case('dryno3')
      allvars(n_vars)%indices(1)=ino3_a
    case('wetno3')
      allvars(n_vars)%indices(1)=ino3_a
    case('dryhno3')
      allvars(n_vars)%indices(1)=ihno3
    case('wethno3')
      allvars(n_vars)%indices(1)=ihno3
    case('dryno2')
      allvars(n_vars)%indices(1)=ino2
    case('dryno')
      allvars(n_vars)%indices(1)=ino
    case('drypan')
      allvars(n_vars)%indices(1)=ipan
    case('sfo3max')
      allvars(n_vars)%indices(1)=io3
    case('mono')
      allvars(n_vars)%indices(1)=iterp
    case('co2mass')
      allvars(n_vars)%indices(1)=-1! not available
    case('chepsoa','chepsoa3d','chepsoa2d')      
      allvars(n_vars)%indices(1)=-1
    case('flashrate')
      allvars(n_vars)%indices(1)=-1
    case('md')
      allvars(n_vars)%indices(1)=-1
    case('sedustCI')      
      allvars(n_vars)%indices(1)=iducoi
    case('depdust')
      allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/)
    case('concdust')
      allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/)
    case('conccn')
      allvars(n_vars)%indices(1:7)=(/inus_n,iais_n,iacs_n,icos_n,iaii_n,iaci_n,icoi_n/)
    case('sconcss')
      allvars(n_vars)%indices(1:2)=(/issacs,isscos/)
    case('loaddust')
      allvars(n_vars)%indices(1:4)=(/iduacs,iducos,iduaci,iducoi/)
    case('lossch4')
      allvars(n_vars)%indices(1)=ich4
    case('lossco')
      allvars(n_vars)%indices(1)=ico

    end select


    if (class=='crescendo') then
      select case (trim(shortname))
      case('od440dust')
        allvars(n_vars)%indices(1)=-1
        if (freq=='day' .and. data_dims==2)then
          od440dustday=n_vars
        end if
      case('od440aer')
        allvars(n_vars)%indices(1)=-1
        if (freq=='1hr') then
          od440aerhr=n_vars
        else if (freq=='day') then
          od440aerdaily=n_vars
        else if (freq=='mon') then
          od440aermonthly=n_vars
        end if
      case('od550aer')
        allvars(n_vars)%indices(1)=-1
        if (freq=='1hr')then
          od550aerhr=n_vars
        end if
      case('od550dust')
        if (freq=='day' .and. data_dims==2)then
          od550dust2dday=n_vars
        end if
      case('od5503ddust')
        if (freq=='day' .and. data_dims==3)then
          od550dust3dday=n_vars
        end if
      case('ec550aer')
        allvars(n_vars)%indices(1)=-1
        ec550aermon=n_vars
      end select

    end if
    if (class=='sf1d') then
      select case (trim(shortname))
      case('od550aer')
        allvars(n_vars)%indices(1)=-1
        od550aerdaily=n_vars
      case('sfo3max')
        allvars(n_vars)%indices(1)=io3
      end select
    end if
    if(class=='optics') then
      select case (trim(shortname))
      case('ec550aer')
        allvars(n_vars)%indices(1)=-1

        ec550aer=n_vars
      case('od550aer')
        allvars(n_vars)%indices(1)=-1
        od550aer=n_vars
      case('od550so4')
        allvars(n_vars)%indices(1)=-1
        od550so4=n_vars
      case('od550bc')
        allvars(n_vars)%indices(1)=-1
        od550bc=n_vars
      case('od550oa')
        allvars(n_vars)%indices(1)=-1
        od550oa=n_vars
      case('od550soa')
        allvars(n_vars)%indices(1)=-1
        od550soa=n_vars
      case('od550ss')
        allvars(n_vars)%indices(1)=-1
        od550ss=n_vars
      case('od550dust')
        allvars(n_vars)%indices(1)=-1
        od550dust=n_vars
      case('od550no3')
        allvars(n_vars)%indices(1)=-1
        od550no3=n_vars
      case('od550aerh2o')
        allvars(n_vars)%indices(1)=-1
        od550aerh2o=n_vars
      case('od550lt1aer')
        allvars(n_vars)%indices(1)=-1
        od550lt1aer=n_vars
      case('abs550aer')
        allvars(n_vars)%indices(1)=-1
        abs550aer=n_vars
      case('od440aer')
        allvars(n_vars)%indices(1)=-1
        if  (freq=='mon') then
          od440aer=n_vars
          !else if (freq=='day') then
          !   od440aerdaily=n_vars
        end if
      case('od350aer')
        allvars(n_vars)%indices(1)=-1
        od350aer=n_vars
      case('od870aer')
        allvars(n_vars)%indices(1)=-1
        od870aer=n_vars
      end select
    end if
    call goLabel()
    status = 0
    call GO_Timer_End( itim_addvar, status )
    IF_NOTOK_RETURN(status=1)
  end subroutine add_variable

  subroutine write_attributes(status,j_var)
    implicit none
    integer     :: j_var
    integer,intent(out)::status
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_writeattr'
    call goLabel(rname)
    call GO_Timer_Start( itim_attr, status )
    IF_NOTOK_RETURN(status=1)
    call MDF_Put_Att( allvars(j_var)%fileunit, MDF_GLOBAL, 'title',               'Model output for Aerchemmip', status )
    IF_NOTOK_MDF(fid= allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'units',          'degrees_east', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'axis',          'X', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'long_name',          'longitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lon_varid , 'standard_name',          'longitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'units',          'degrees_north', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'axis',          'Y', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'long_name',          'latitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lat_varid , 'standard_name',          'latitude', status)
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    ! allvars(j_var)%lev_varid
    if (allvars(j_var)%dims==3) then
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'units',          'level', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'axis',          'Z', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'long_name',          'hybrid model level at layer midpoints', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'standard_name',          'hybrid_model_level', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'formula',          'a+b*ps', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%lev_varid , 'positive',          'up', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    


      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyam_varid , 'long_name',          'hybrid A coefficient at layer midpoints', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyam_varid , 'units',          'Pa', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybm_varid , 'long_name',          'hybrid B coefficient at layer midpoints', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybm_varid , 'units',          '1', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyai_varid , 'long_name',          'hybrid A coefficient at layer interfaces', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hyai_varid , 'units',          'Pa', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybi_varid , 'long_name',          'hybrid B coefficient at layer interfaces', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%hybi_varid , 'units',          '1', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    END if
    if (allvars(j_var)%table_id/='AERfx')then
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'units',          'days since 1750-01-01 00:00:00', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      !    call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'axis',          'X', status)
      !    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'calendar',          'gregorian', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'long_name',          'time', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'axis',          'T', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    

      !time bounds
      call MDF_Put_Att( allvars(j_var)%fileunit,allvars(j_var)%time_varid , 'bounds',          'time_bounds', status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)    
    end if

    !experiment=
    !CMIP6 global attributes:
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,    'Conventions', 'CF-1.7 CMIP-6.0 UGRID-0.9', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'activity_id',  'AerChemMIP', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'branch_method','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'creation_date','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'data_specs_version','1.0.0', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'experiment',trim(experiment), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'experiment_id',trim(experiment_id), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'forcing_index',trim(forcing_i), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'frequency',trim(allvars(j_var)%freq), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'further_info_url','MISSING', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'grid','native '//trim(grid)//' degree grid', status)!module variables
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'grid_label',trim(grid_label), status)!module variables
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'initialization_index',trim(initialization_i), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'institution','KNMI, The Netherlands; SMHI, Sweden; DMI, Denmark; &
         &AEMET, Spain; Met Eireann, Ireland; CNR-ISAC, Italy; Instituto de Meteorologia, Portugal; FMI, Finland; BSC, Spain; &
         &Centro de Geofisica, University of Lisbon, Portugal; ENEA, Italy; Geomar, Germany; Geophysical Institute, University of Bergen, Norway; &
         &ICHEC, Ireland; ICTP, Italy; IMAU, The Netherlands; IRV, Sweden;  Lund University, Sweden; &
         &Meteorologiska Institutionen, Stockholms University, Sweden; Niels Bohr Institute, University of Copenhagen, Denmark; &
         &NTNU, Norway; SARA, The Netherlands; Unite ASTR, Belgium; Universiteit Utrecht, The Netherlands; &
         &Universiteit Wageningen, The Netherlands; University College Dublin, Ireland; Vrije Universiteit Amsterdam, the Netherlands; &
         &University of Helsinki, Finland; KIT, Karlsruhe, Germany; USC, University of Santiago de Compostela, Spain; &
         &Uppsala Universitet, Sweden; NLeSC, Netherlands eScience Center, The Netherlands', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'institution_id','EC-Earth-Consortium', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'license','NEEDS DEFINING', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'mip_era','CMIP6', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'nominal_resolution','250 km', status)!dmax
    !dmax=r*phi/2*(1+((phi**2+lamb**2)/(phi*lamb))*np.arctan(lamb/phi))=348 r=6371, phi=2(lat), lamb=3(long)
    !CMIP6 global attributes: 100 < dmax < 350 -> nominal resolution 250 km
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'physics_index',trim(physics_i), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'product','output', status)!only choice
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'realization_index',trim(realization_i), status)!1 for primary or single realization
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'realm',trim(realm), status)!  depends on run, some are AGCM
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'source',trim(source), status)!
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'source_id',trim(source_id), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'source_type',trim(source_type), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'sub_experiment','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'sub_experiment_id','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'table_id',trim(allvars(j_var)%table_id), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'tracking_id','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'variable_id',trim(allvars(j_var)%vname), status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'variant_label','', status)
    call MDF_Put_Att( allvars(j_var)%fileunit,MDF_GLOBAL ,        'variant_label','', status)

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

    call goLabel() 
    status = 0

  end subroutine write_attributes

  subroutine write_dimensions(status,j_var)
    use dims,           only : iyear0          !current year
    use go_date,        only : days_in_year,newDate
    use partools    , only : isRoot,myid
    implicit none
    integer     :: j_var
    integer,intent(out)::status
    integer     :: i1,i2,j1,j2,imr,jmr,lmr
    integer     :: lon_varid,lonid,lon_dimid
    integer     :: lat_varid,latid,lat_dimid
    integer     :: lev_varid,levid,lev_dimid
    integer     :: hym_varid,hym_dimid
    integer     :: hyi_varid,hyi_dimid
    integer     :: time_varid,timeid,time_dimid,bounds_dimid,bounds_varid,boudid
    ! most of data  is monthly mean, but change to dynamic number of output steps needed
    integer     :: nout_steps!=12 
    integer     :: nhym
    integer     :: nhyi 

    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_write_dim'
    call goLabel(rname)

    ! define dimensions
    !call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=dimension_data%nlon
    jmr=dimension_data%nlat
    lmr=dimension_data%nlev
    nhym=lmr
    nhyi=lmr+1

    ! With parallel netcdf whole netcdf must be reserved at the time of initialisation
    ! therefore we need to know the number of output steps per file.

    ! Define number of output steps in a file depending on the output frequency
    ! use newdate to create TDate structure, and use that in days_in_year
    if (allvars(j_var)%table_id=='AERhr')then
      nout_steps=24*days_in_year(newdate(iyear0))
    else if (allvars(j_var)%table_id=='AER6hr')then
      nout_steps=4*days_in_year(newdate(iyear0))
    else if (allvars(j_var)%table_id=='AERday')then
      nout_steps=days_in_year(newdate(iyear0))
    else if (allvars(j_var)%table_id=='AERmon'.or. (allvars(j_var)%table_id=='AERmonZ'))then
      nout_steps=12
    end if

    if (isroot) then
      !DEFINE DIMENSIONS
      call MDF_Def_Dim( allvars(j_var)%fileunit, 'lon',     imr,lon_dimid, status )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      call MDF_Def_Dim( allvars(j_var)%fileunit, 'lat',      jmr, lat_dimid, status )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      if (allvars(j_var)%dims==3) then
        if (trim(allvars(j_var)%vname)=='phalf') then
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'lev',        lmr+1, lev_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        else 
          call MDF_Def_Dim( allvars(j_var)%fileunit, 'lev',        lmr, lev_dimid, status )
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        end if
      end if

      if (allvars(j_var)%table_id/='AERfx')then
        !call MDF_Def_Dim( allvars(j_var)%fileunit, 'time',       nout_steps, time_dimid, status )
        call MDF_Def_Dim( allvars(j_var)%fileunit, 'time',       MDF_UNLIMITED, time_dimid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        call MDF_Def_Dim( allvars(j_var)%fileunit, 'bounds',       2, bounds_dimid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      end if
      if (allvars(j_var)%dims==3) then
        call MDF_Def_Dim( allvars(j_var)%fileunit, 'nhym',       nhym, hym_dimid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        call MDF_Def_Dim( allvars(j_var)%fileunit, 'nhyi',       nhyi, hyi_dimid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      end if

      ! define dimension variables 
      ! longitude
      call MDF_Def_Var( allvars(j_var)%fileunit, 'lon', MDF_DOUBLE, &
           (/ lon_dimid/), allvars(j_var)%lon_varid, status )


      ! define latitude variable
      call MDF_Def_Var( allvars(j_var)%fileunit, 'lat', MDF_DOUBLE, &
           (/ lat_dimid/), allvars(j_var)%lat_varid, status )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
      ! level
      if (allvars(j_var)%dims==3) then
        call MDF_Def_Var( allvars(j_var)%fileunit, 'lev', MDF_DOUBLE, &
             (/ lev_dimid/), allvars(j_var)%lev_varid, status )
      end if

      if (allvars(j_var)%table_id/='AERfx')then

        call MDF_Def_Var( allvars(j_var)%fileunit, 'time', MDF_DOUBLE, &
             (/ time_dimid/), allvars(j_var)%time_varid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
        call MDF_Def_Var( allvars(j_var)%fileunit, 'time_bounds', MDF_DOUBLE, &
             (/ bounds_dimid,time_dimid/), allvars(j_var)%bounds_varid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
      end if
      if (allvars(j_var)%dims==3) then
        ! define hybm variable
        call MDF_Def_Var( allvars(j_var)%fileunit, 'hybm', MDF_DOUBLE, &
             (/ hym_dimid/), allvars(j_var)%hybm_varid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
        ! define hyam variable
        call MDF_Def_Var( allvars(j_var)%fileunit, 'hyam', MDF_DOUBLE, &
             (/ hym_dimid/), allvars(j_var)%hyam_varid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
        ! define hybi variable
        call MDF_Def_Var( allvars(j_var)%fileunit, 'hybi', MDF_DOUBLE, &
             (/ hyi_dimid/), allvars(j_var)%hybi_varid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
        ! define hyai variable
        call MDF_Def_Var( allvars(j_var)%fileunit, 'hyai', MDF_DOUBLE, &
             (/ hyi_dimid/), allvars(j_var)%hyai_varid, status )
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit) 
      end if

      ! * Write out dimension variable values

      ! Write out hybm 
      if (allvars(j_var)%dims==3) then
        ! midpoints
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hybm_varid,levi%fb, status)
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        ! Write out hyam
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hyam_varid,levi%fa, status)
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        ! interfaces 
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hybi_varid,levi%b, status)
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        ! Write out hyai
        if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%hyai_varid,levi%a, status)
        IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
      end if

      ! Write out the longitudes
      if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lon_varid, dimension_data%lon, status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

      !write out the latitude to variable
      if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit,  allvars(j_var)%lat_varid, dimension_data%lat, status)
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

      if (allvars(j_var)%dims==3) then
        if (trim(allvars(j_var)%vname)=='phalf') then
          ! Write out the levels
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/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,35/), status)
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
        else
          if (isroot)call MDF_Put_Var( allvars(j_var)%fileunit, allvars(j_var)%lev_varid, (/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)
          IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

        end if

      end if
      ! Time will be written during output write -steps

    end if
    call goLabel()
    status = 0

  end subroutine write_dimensions

  subroutine write_var(status,j_var)
    implicit none
    integer     :: j_var
    integer,intent(out)::status
    integer     :: i1,i2,j1,j2,imr,jmr,lmr
    integer     :: lon_varid,lonid
    integer     :: lat_varid,latid
    integer     :: lev_varid,levid
    integer     :: time_varid,timeid


    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_write_dim'
    call goLabel(rname)

    ! define dimensions
    !call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    imr=dimension_data%nlon
    jmr=dimension_data%nlat
    lmr=dimension_data%nlev
    ! define dimension variables 
    ! longitude

    if (allvars(j_var)%dims==2.and.allvars(j_var)%table_id/='AERfx') then 
      call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, &
           (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid, allvars(j_var)%time_varid/), allvars(j_var)%varid, status )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

    else  if (allvars(j_var)%dims==2.and.allvars(j_var)%table_id=='AERfx') then 
      call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, &
           (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid/), allvars(j_var)%varid, status )
      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    else
      !(allvars(j_var)%dims==3)


      call MDF_Def_Var( allvars(j_var)%fileunit, allvars(j_var)%vname, MDF_DOUBLE, &
           (/allvars(j_var)%lon_varid, allvars(j_var)%lat_varid, allvars(j_var)%lev_varid,allvars(j_var)%time_varid/), allvars(j_var)%varid, status )

      IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    end if
    ! Write out the longitudes


    call MDF_Put_Att( allvars(j_var)%fileunit, allvars(j_var)%varid, 'long_name',     trim(allvars(j_var)%lname), status )
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att(allvars(j_var)%fileunit,allvars(j_var)%varid, 'standard_name', trim(allvars(j_var)%standard_name), status )
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)
    call MDF_Put_Att(allvars(j_var)%fileunit , allvars(j_var)%varid, 'units',         trim(allvars(j_var)%unit), status )
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

    call MDF_EndDef( allvars(j_var)%fileunit, status )
    IF_NOTOK_MDF(fid=allvars(j_var)%fileunit)

    call goLabel()
    status = 0

  end subroutine write_var

  subroutine calculate_optics( region, dhour, l_do_ec550aer_only,status)
    use optics,      only : optics_aop_get
    use MeteoData,   only : gph_dat,temper_dat
    USE toolbox,       only : ltropo_ifs,     lvlpress
    use dims,        only : lm  
    implicit none

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


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

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

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

    integer :: ncontr, lvec, lct, l, indoffset, nwl
    integer :: nadd, nadd_max, indoffset_stat, indabs
    integer :: iadd
    integer :: region,dhour,status
    real    :: no3fac, wght, dwght, tabs
    real,dimension(:),allocatable      :: tx
    integer                            :: itropo
    real, dimension(:,:,:), pointer    :: gph   ! height (incl. oro)
    real,dimension(:,:,:),  pointer    :: t           ! temperature (K)
    logical :: l_do_ec550aer_only
    character(len=*), parameter ::  rname = mname//'/output_aerchemmip_optics'
    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

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

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

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

    ! Total, SO4, BC, OC, SS, DU, NO3, Water, Fine, Fine Dust, Fine SS

    ! Total(water), Total(nowater)
    !ncontr = 10
    ncontr = 12
    ! Total, SO4, BC, OC,SOA, SS, DU, NO3, Water, lt1aer

    dtime=dhour*3600

    !TB Following additional fields are nnecessary but call to optics needs it...

    ! 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, .false., 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/) )
    end do ! lwl

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


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

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

    ! cycle over wavelengths
    do lwl = 1, size(wdep_out) 
      if (.not. l_do_ec550aer_only)then 
        ! if split and if wavelength=550nm
        if( wdep_out(lwl)%split ) then 
          if (wdep_out(lwl)%wl == 0.550) then 
            ! for 550nm: 
            ! 1) AODs 
            !    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))
            allocate(tempdust2dday(i1:i2,j1:j2))

            do lct = 1, ncontr
              temp2d = 0.0
              tempdust2dday=0.0
              do j = j1, j2
                do i = i1, i2
                  ! multiply with height elements and sum up
                  tabs = 0.0
                  ! AOD output will only be for troposphere in data request
                  tx(:)=t(i,j,:)
                  itropo = ltropo_ifs(region,i,j,tx,lm(region))
                  do l = 1, itropo!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 ! TOTAL AOD
                      ! 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))


                    else if (lct==4) then ! OAAOD
                      ! add SOA aod (5)to POM aod (4) (AerChemMIP expects total oa in aod550oa)
                      temp2d(i,j) = temp2d(i,j) + opt_ext(i,j,l,lwl,5) * & 
                           (gph_dat(region)%data(i,j,l+1) - & 
                           gph_dat(region)%data(i,j,l  ))
                    else if (lct==7) then ! DUSTAOD
                      if ( wdep_out(lwl)%wl == 0.550) then
                        if (crescendo_out) allvars(od550dust3dday)%data3D(i,j,l)= opt_ext(i,j,l,lwl,lct) * & 
                             (gph_dat(region)%data(i,j,l+1) - & 
                             gph_dat(region)%data(i,j,l  ))

                        tempdust2dday(i,j)=tempdust2dday(i,j)+ opt_ext(i,j,l,lwl,lct) * & 
                             (gph_dat(region)%data(i,j,l+1) - & 
                             gph_dat(region)%data(i,j,l  ))
                      else 
                        tempdust2dday(i,j)=tempdust2dday(i,j)+ opt_ext(i,j,l,lwl,lct) * & 
                             (gph_dat(region)%data(i,j,l+1) - & 
                             gph_dat(region)%data(i,j,l  ))
                      end if
                    end if
                  end do
                  ! Absorption: do this only once for the totals
                  if( lct == 1 ) then 

                    allvars(abs550aer)%data2D(i,j)=allvars(abs550aer)%data2D(i,j) + tabs*dtime
                    !extra output for total od550aer (ncontr==1)
                    allvars(od550aerdaily)%data2D(i,j)=allvars(od550aerdaily)%data2D(i,j)  + temp2d(i,j)*dtime

                    if (crescendo_out)then
                      allvars(od550aerhr)%data2D(i,j)= dtime*temp2d(i,j)
                    end if
                    allvars(od550aer)%data2D(i,j)=allvars(od550aer)%data2D(i,j)  + temp2d(i,j)*dtime
                  else if (lct<11)Then !AOD by compounds 

                    allvars(indoffset+lct-1)%data2D(i,j)=allvars(indoffset+lct-1)%data2D(i,j) + temp2d(i,j)*dtime

                    if (crescendo_out .and. lct==7 .and. wdep_out(lwl)%wl == 0.550) then !DUSTAOD for 550nm
                      allvars(od550dust2dday)%data2D(i,j)=allvars(od550dust2dday)%data2D(i,j)  + tempdust2dday(i,j)*dtime
                    end if

                  end if

                end do
              end do

            end do
            deallocate( temp2d )
            deallocate( tempdust2dday)
            !if 440 has splits do this
          else  if (wdep_out(lwl)%wl == 0.440 ) then
            indoffset = od440aer
            !not needed for AERCHEMMIP
            indabs = -1
            !abs440aer
            ! fill AODs
            allocate(tempdust2d440day(i1:i2,j1:j2))
            allocate(temp2d(i1:i2,j1:j2))
            tempdust2d440day=0.0
            temp2d = 0.0
            tempdust2d440day=0.0
            do j = j1, j2
              do i = i1, i2
                ! AOD output will only be for troposphere in data request
                tx(:)=t(i,j,:)
                itropo = ltropo_ifs(region,i,j,tx,lm(region))

                ! multiply with height elements and sum up
                tabs = 0.0

                do l = 1, itropo!lmr
                  ! od440aer
                  lct=1 ! total aod

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

                  !OD440DUST

                  lct=7
                  tempdust2d440day(i,j)=tempdust2d440day(i,j)+ opt_ext(i,j,l,lwl,lct) * & 
                       (gph_dat(region)%data(i,j,l+1) - & 
                       gph_dat(region)%data(i,j,l  ))  

                end do
                if (indabs>0) then 
                  allvars(indabs)%data2D(i,j)=allvars(indabs)%data2D(i,j) + tabs*dtime
                end if
              end do
            end do

            allvars(od440aer)%data2D=allvars(od440aer)%data2D + temp2d*dtime
            if (crescendo_out)then
              allvars(od440dustday)%data2D=allvars(od440dustday)%data2D + tempdust2d440day*dtime
              allvars(od440aerdaily)%data2D=allvars(od440aerdaily)%data2D + temp2d*dtime
              allvars(od440aerhr)%data2D=temp2d*dtime
            end if

            deallocate( tempdust2d440day)
            deallocate( temp2d )
          end if





        else !NON SPLITS
          ! 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.870 ) then
            indoffset = od870aer
            !not needed for AERCHEMMIP
            indabs = -1
            !abs870aer
          elseif ( wdep_out(lwl)%wl == 0.440 ) then
            indoffset = -1 !od440aer
            !not needed for AERCHEMMIP
            indabs = -1
            !abs350aer
          elseif ( wdep_out(lwl)%wl == 0.350 ) then
            indoffset = od350aer
            !not needed for AERCHEMMIP
            indabs = -1
            !abs350aer
          else
            write(gol,*) 'user_output_aerchemmip: wrong wavelength given for aerocom diagnostics' ; call goErr
            TRACEBACK
            status=1; return
          end if

          ! fill AODs, total only
          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
              if (indabs>0) then 
                allvars(indabs)%data2D(i,j)=allvars(indabs)%data2D(i,j) + tabs*dtime
              end if

            end do
          end do

          if ((indoffset)==od870aer) then
            allvars(od870aer)%data2D=allvars(od870aer)%data2D + temp2d*dtime
          end if

          deallocate(temp2d)

        endif
      end if
      ! Extinction in 3D
      if ( wdep_out(lwl)%wl == 0.550 ) then
        allvars(ec550aer)%data3D= opt_ext(:,:,:,lwl,1)
        if (crescendo_out) allvars(ec550aermon)%data3D=allvars(ec550aermon)%data3D + opt_ext(:,:,:,lwl,1)*dtime

      endif


    end do

    deallocate( opt_ext, opt_a, opt_g, opt_add )
    deallocate( tx )

    call goLabel() ; status = 0

  end subroutine calculate_optics



  real function mode_fraction(r,limit,imode) result(modfrac)
    Use mo_aero_m7,  only : nmod, cmedr2mmedr, sigma
    implicit none
    !real :: var
    real :: r
    real :: z
    real :: limit
    ! real :: sigma
    real :: hr2=0.5*sqrt(2.0)
    integer::imode
    z=0.0
    if( r<100*tiny(1.0))then
      modfrac=1.0
    else
      z=(log(limit)-log(r*cmedr2mmedr(imode)))/log(sigma(imode))

      modfrac=0.5+0.5*erf(z*hr2)
    end if
  end function mode_fraction
end MODULE USER_OUTPUT_AERCHEMMIP