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
  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
  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), 'mol s-1', 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.

  !SPECIAL
  !HOURLY
  character(len=8),dimension(1),parameter:: ps6hr=(/'ps'/)
  character(len=11),dimension(1),parameter:: ps6hrunit=(/units(ipresunit)/)

  character(len=8),dimension(6),parameter:: hourly_var=(/'ps', 'sfno2', 'sfo3', 'sfpm25', 'tas', 'ec550aer'/)
  character(len=11),dimension(6),parameter:: hourly_varunit=(/units(ipresunit), units(ivmrunit), units(ivmrunit), &
       units(immrunit), units(itempunit), units(iextunit)/)
  !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 :: nlev_region ! also global
    !integer :: nlev_region ! also global
    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 usest 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))
       !imr = i2-i1+1
       !jmr = j2-j1+1       
       ! 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
    
    !Initialise 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(ps6hr)
       call add_variable(-1,trim(ps6hr(i)),trim(ps6hr(i)),ps6hrunit(i),2,status,'ps6h','AER6hr',-1.0)
    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
    do i=1,size(hourly_var)
              
       if (trim(hourly_var(i))=='ec550aer')then
          call add_variable(-1,trim(hourly_var(i)),'optics '//trim(hourly_var(i)), units(iextunit),3,status,'optics','AER6hr',-1.0)
       else
          call add_variable(-1,trim(hourly_var(i)),trim(hourly_var(i)),hourly_varunit(i),2,status,'sf1h','AERhr',-1.0)
       end if
    end do
    ! 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

    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)
    ! 
   
    do j_var = 1, n_vars
       ! initialise a single file for each variables as per AERCHEMMIP specification
       
       ! 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)
       
       
       !For each file 
       ! 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)
       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
             ! Gridbox area
             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
             ! Gridbox area

             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')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'
    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
                   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
                   if (tracers(i,j,1,io3)>  allvars(j_var)%data2D(i,j)) then
                      allvars(j_var)%data2D(i,j)=tracers(i,j,1,io3)*xmair/xmo3/m_dat(region)%data(i,j,1) !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