Index: /branches/flexpart91_hasod/src_parallel/FLEXPART.f90
===================================================================
--- /branches/flexpart91_hasod/src_parallel/FLEXPART.f90	(revision 8)
+++ /branches/flexpart91_hasod/src_parallel/FLEXPART.f90	(revision 9)
@@ -50,5 +50,5 @@
 
 #ifdef NETCDF_OUTPUT
-  use ncdf_output_mod, only: writeheader_ncdf
+  use nc_output_mod, only: writeheader_ncdf
 #endif
 
Index: /branches/flexpart91_hasod/src_parallel/Makefile.netcdf
===================================================================
--- /branches/flexpart91_hasod/src_parallel/Makefile.netcdf	(revision 9)
+++ /branches/flexpart91_hasod/src_parallel/Makefile.netcdf	(revision 9)
@@ -0,0 +1,97 @@
+SHELL = /bin/bash
+MAIN = FLEXPART_GFORTRAN
+#
+
+FC       = gfortran
+#INCPATH  = /usr/local/ecmwf_tools/reloc/include
+#LIBPATH1 = /usr/local/ecmwf_tools/reloc/lib
+INCPATH  = /home/sharald/local/include
+NCINCPATH = /home/sharald/local/include
+LIBPATH1 = /home/sharald/local/lib64
+LIBPATH2 = /usr/local/lib
+NCLIBPATH = /home/sharald/local/lib64
+#FFLAGS   = -O3 -g -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+#FFLAGS   = -O0 -g -cpp -fopenmp -m64 -fbacktrace -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+FFLAGS   = -O2 -cpp -fopenmp -m64 -fstack-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -DNETCDF_OUTPUT -I$(INCPATH) -I$(NCINCPATH)
+#FFLAGS   = -O1 -g -cpp -m64 -fbacktrace -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+LDFLAGS  = $(FFLAGS) -L$(LIBPATH1) -L$(LIBPATH2) -L$(NCLIBPATH) -lnetcdf -lnetcdff -lgrib_api_f90 -lgrib_api -lm -lopenjpeg
+#
+
+MODOBJS = \
+par_mod.o               com_mod.o \
+conv_mod.o              hanna_mod.o \
+interpol_mod.o          cmapf_mod.o \
+unc_mod.o               oh_mod.o \
+xmass_mod.o             flux_mod.o \
+point_mod.o             outg_mod.o \
+random_mod.o 		nc_output_mod.o
+
+OBJECTS = \
+writeheader.o           assignland.o\
+calcpar.o               part0.o \
+caldate.o               partdep.o \
+coordtrafo.o            psih.o \
+raerod.o \
+drydepokernel.o          \
+readavailable.o \
+ew.o                    readcommand.o \
+advance.o               readdepo.o \
+releaseparticles.o      psim.o     \
+FLEXPART.o              readlanduse.o \
+getfields.o             init_domainfill.o\
+interpol_wind.o         readoutgrid.o \
+interpol_all.o          readpaths.o \
+getrb.o                 readreceptors.o \
+getrc.o                 readreleases.o \
+getvdep.o               readspecies.o \
+interpol_misslev.o      readwind.o \
+conccalc.o              richardson.o \
+concoutput.o            scalev.o \
+pbl_profile.o           readOHfield.o\
+juldate.o               timemanager.o \
+interpol_vdep.o         interpol_rain.o \
+verttransform.o         partoutput.o \
+hanna.o                 wetdepokernel.o \
+mean.o                  wetdepo.o \
+hanna_short.o           windalign.o \
+obukhov.o               gridcheck.o \
+hanna1.o                initialize.o \
+                        gridcheck_nests.o \
+readwind_nests.o        calcpar_nests.o \
+verttransform_nests.o   interpol_all_nests.o \
+interpol_wind_nests.o   interpol_misslev_nests.o \
+interpol_vdep_nests.o   interpol_rain_nests.o \
+getvdep_nests.o \
+readageclasses.o        readpartpositions.o \
+calcfluxes.o            fluxoutput.o \
+qvsat.o                 skplin.o \
+convmix.o               calcmatrix.o \
+convect43c.o               redist.o \
+sort2.o                 distance.o \
+centerofmass.o          plumetraj.o \
+openouttraj.o           calcpv.o \
+calcpv_nests.o          distance2.o \
+clustering.o            interpol_wind_short.o \
+interpol_wind_short_nests.o shift_field_0.o \
+shift_field.o           outgrid_init.o \
+openreceptors.o         boundcond_domainfill.o\
+partoutput_short.o      readoutgrid_nest.o \
+outgrid_init_nest.o     writeheader_nest.o \
+concoutput_nest.o       wetdepokernel_nest.o \
+drydepokernel_nest.o    zenithangle.o \
+ohreaction.o            getvdep_nests.o \
+initial_cond_calc.o     initial_cond_output.o \
+dynamic_viscosity.o     get_settling.o
+
+
+$(MAIN): $(MODOBJS) $(OBJECTS)
+	$(FC) *.o -o $(MAIN) $(LDFLAGS)
+
+$(OBJECTS): $(MODOBJS)
+
+%.o: %.f90
+	$(FC) -c $(FFLAGS) $<
+
+clean:
+	rm *.o *.mod
+
Index: /branches/flexpart91_hasod/src_parallel/nc_output_mod.f90
===================================================================
--- /branches/flexpart91_hasod/src_parallel/nc_output_mod.f90	(revision 9)
+++ /branches/flexpart91_hasod/src_parallel/nc_output_mod.f90	(revision 9)
@@ -0,0 +1,1295 @@
+!**********************************************************************
+! Copyright 2013                                                      *
+! Dominik Brunner                                                     *
+!                                                                     *
+! This file is part of FLEXPART-COSMO                                 *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+  !*****************************************************************************
+  !                                                                            *
+  !  This module handles all gridded netcdf output for concentration or        *
+  !  residence time and wet and dry deposition output.                         *
+  !                                                                            *
+  !  - writeheader_ncdf generates file including all information previously    *
+  !    stored in separate header files                                         *
+  !  - concoutput_ncdf write concentration output and wet and dry deposition   *
+  !                                                                            *
+  !     Author: D. Brunner                                                     *
+  !                                                                            *
+  !     12 April 2013                                                          *
+  !                                                                            *
+  !*****************************************************************************
+
+#ifdef NETCDF_OUTPUT
+
+module nc_output_mod
+
+  use netcdf
+
+  use point_mod, only: ireleasestart,ireleaseend,kindz,&
+                       xpoint1,ypoint1,xpoint2,ypoint2,zpoint1,zpoint2,npart,xmass
+  use outg_mod,  only: outheight,oroout,densityoutgrid,factor3d,volume,&
+                       wetgrid,wetgridsigma,drygrid,drygridsigma,grid,gridsigma,&
+                       area,arean,volumen, orooutn
+  use par_mod,   only: dp, maxspec, maxreceptor, nclassunc,&
+                       unitoutrecept,unitoutreceptppt, nxmax
+  use com_mod,   only: path,length,ldirect,ibdate,ibtime,iedate,ietime, &
+                       loutstep,loutaver,loutsample,outlon0,outlat0,&
+                       numxgrid,numygrid,dxout,dyout,numzgrid, height, &
+                       outlon0n,outlat0n,dxoutn,dyoutn,numxgridn,numygridn, &
+                       nspec,maxpointspec_act,species,numpoint,&
+                       dx,xlon0,dy,ylat0,compoint,method,lsubgrid,lconvection,&
+                       ind_source,ind_receptor,nageclass,lage,&
+                       drydep,wetdep,decay,weta,wetb, numbnests, &
+                       reldiff,henry,f0,density,dquer,dsigma,dryvel,&
+                       weightmolar,ohreact,spec_ass,kao,vsetaver,&
+                       ! for concoutput_ncdf and concoutput_nest_ncdf
+                       nxmin1,nymin1,nz,oro,oron,rho,rhon,&
+                       memind,xresoln,yresoln,xrn, xln, yrn,yln,nxn,nyn,&
+                       xreceptor,yreceptor,numreceptor,creceptor,iout, &
+                       itsplit, lsynctime, ctl, ifine, lagespectra, ipin, &
+                       ioutputforeachrelease, iflux, mdomainfill, mquasilag, & 
+                       nested_output, ipout
+
+  implicit none
+
+  ! parameter for data compression (1-9, 9 = most aggressive)
+  integer, parameter :: deflate_level = 9
+  logical, parameter :: min_size = .false.   ! if set true, redundant fields (topography) are not written to minimize file size
+  character(len=255), parameter :: institution = 'Empa, Duebendorf, Switzerland'
+
+  integer            :: tpointer
+  character(len=255) :: ncfname, ncfnamen
+
+  ! netcdf dimension and variable IDs for main and nested output grid
+  integer, dimension(maxspec) :: specID,specIDppt, wdspecID,ddspecID
+  integer, dimension(maxspec) :: specIDn,specIDnppt, wdspecIDn,ddspecIDn
+  integer                     :: timeID, timeIDn
+  integer, dimension(6)       :: dimids, dimidsn
+  integer, dimension(5)       :: depdimids, depdimidsn
+  real,parameter :: eps=nxmax/3.e5
+
+  private:: writemetadata, output_units, nf90_err
+
+contains
+
+!****************************************************************
+! determine output units (see table 1 in Stohl et al., ACP 2005
+!****************************************************************
+subroutine output_units(units)
+  character(len=15), intent(out) :: units
+  if (ldirect.eq.1) then
+     ! forward simulation
+     if (ind_source.eq.1) then
+        if (ind_receptor.eq.1) then
+           units = 'ng m-3'   ! hes the kg in Tab1 is only indicating the units of the relase not the output
+        else
+           units = 'ng kg-1'
+        endif
+     else
+        if (ind_receptor.eq.1) then
+           units = 'ng m-3'
+        else
+           units = 'ng kg-1'
+        endif
+     endif
+  else
+     ! backward simulation
+     if (ind_source.eq.1) then
+        if (ind_receptor.eq.1) then
+           units = 's'
+        else 
+           units = 's m3 kg-1'
+        endif
+     else
+        if (ind_receptor.eq.1) then
+           units = 's kg m-3'
+        else
+           units = 's'
+        endif
+     endif
+  endif
+end subroutine output_units
+
+
+!****************************************************************
+! write metadata to netCDF file 
+!****************************************************************
+subroutine writemetadata(ncid,lnest)
+
+  integer, intent(in) :: ncid
+  logical, intent(in) :: lnest
+  integer             :: status
+  character           :: time*10,date*8,adate*8,atime*6
+  character(5)        :: zone
+  character(255)      :: login_name, host_name
+
+! gather system information 
+  call date_and_time(date,time,zone)
+  call getlog(login_name)
+  call hostnm(host_name)
+  
+! hes CF convention requires these attributes
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'Conventions', 'CF-1.6'))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'title', 'FLEXPART model output'))
+! call nf90_err(nf90_put_att(ncid, nf90_global, 'institution', trim(institution)))
+! call nf90_err(nf90_put_att(ncid, nf90_global, 'source', trim(version_str)//' model output'))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'history', date(1:4)//'-'//date(5:6)// &
+       '-'//date(7:8)//' '//time(1:2)//':'//time(3:4)//' '//zone//'  created by '//  &
+       trim(login_name)//' on '//trim(host_name)))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'references', &
+       'Stohl et al., Atmos. Chem. Phys., 2005, doi:10.5194/acp-5-2461-200'))
+
+  ! attributes describing model run
+  !************************************************************************************
+
+  if (lnest) then
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'outlon0', outlon0n))
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'outlat0', outlat0n))
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'dxout', dxoutn))
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'dyout', dyoutn))
+  else
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'outlon0', outlon0))
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'outlat0', outlat0))
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'dxout', dxout))
+     call nf90_err(nf90_put_att(ncid, nf90_global, 'dyout', dyout))
+  endif
+!	vertical levels stored in grid structure
+
+  ! COMMAND file settings
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ldirect', ldirect))
+  write(adate,'(i8.8)') ibdate
+  write(atime,'(i6.6)') ibtime
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ibdate', adate))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ibtime', atime))
+  write(adate,'(i8.8)') iedate
+  write(atime,'(i6.6)') ietime
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'iedate', adate))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ietime', atime))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'loutstep', loutstep))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'loutaver', loutaver))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'loutsample', loutsample))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'itsplit', itsplit))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'lsynctime', lsynctime))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ctl', ctl))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ifine', ifine))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'iout', iout))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ipout', ipout))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'lsubgrid', lsubgrid))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'lconvection', lconvection))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'lagespectra', lagespectra))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ipin', ipin))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ioutputforeachrelease', ioutputforeachrelease))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'iflux', iflux))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'mdomainfill', mdomainfill))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ind_source', ind_source))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'ind_receptor', ind_receptor))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'mquasilag', mquasilag))
+  call nf90_err(nf90_put_att(ncid, nf90_global, 'nested_output', nested_output))
+ 
+end subroutine writemetadata
+
+
+!****************************************************************
+! netcdf error message handling
+!****************************************************************
+subroutine nf90_err(status)
+  integer, intent (in) :: status
+   if(status /= nf90_noerr) then
+      print *, trim(nf90_strerror(status))
+      stop 'Stopped'
+    end if
+end subroutine nf90_err
+
+
+!****************************************************************
+! Create netcdf file and write header/metadata information
+! lnest = .false. : Create main output file
+! lnest = .true.  : Create nested output file
+!****************************************************************
+subroutine writeheader_ncdf(lnest)
+
+  implicit none
+
+  logical, intent(in) :: lnest
+
+  integer :: ncid, sID, wdsID, ddsID
+  integer :: timeDimID, latDimID, lonDimID, levDimID
+  integer :: nspecDimID, npointDimID, nageclassDimID, ncharDimID
+  integer :: tID, lonID, latID, levID, poleID, lageID, oroID
+  integer :: rellng1ID, rellng2ID, rellat1ID, rellat2ID, relzz1ID, relzz2ID
+  integer :: relcomID, relkindzID, relstartID, relendID, relpartID, relxmassID
+  integer :: nnx, nny 
+  integer, dimension(6)       :: dIDs
+  integer, dimension(5)       :: depdIDs
+  character(len=255)          :: fname
+  character(len=15)           :: units
+  character(len=10)           :: fprefix
+  character(len=3)            :: anspec
+  character                   :: adate*8,atime*6,timeunit*32
+  real, dimension(1000)       :: coord
+
+  integer                     :: cache_size
+  integer, dimension(6)       :: chunksizes
+  integer, dimension(5)       :: dep_chunksizes
+
+  integer                     :: i,ix,jy
+
+  !************************
+  ! Create netcdf file
+  !************************
+
+  if (ldirect.eq.1) then
+     write(adate,'(i8.8)') ibdate
+     write(atime,'(i6.6)') ibtime
+     fprefix = 'grid_conc_'
+  else
+     write(adate,'(i8.8)') iedate
+     write(atime,'(i6.6)') ietime
+     fprefix = 'grid_time_'
+  endif
+
+  if (lnest) then
+     fname = path(2)(1:length(2))//fprefix//adate//atime//'_nest.nc'
+     ncfnamen = fname
+     nnx = numxgridn
+     nny = numygridn
+  else
+     fname = path(2)(1:length(2))//fprefix//adate//atime//'.nc'
+     ncfname = fname
+     nnx = numxgrid
+     nny = numygrid
+  endif
+
+  cache_size = 16 * nnx * nny * numzgrid
+
+  ! setting cache size in bytes. It is set to 4 times the largest data block that is written
+  !   size_type x nx x ny x nz
+  ! create file
+
+  call nf90_err(nf90_create(trim(fname), cmode = nf90_hdf5, ncid = ncid, &
+    cache_size = cache_size))  
+
+  ! create dimensions:
+  !*************************
+  ! time
+  call nf90_err(nf90_def_dim(ncid, 'time', nf90_unlimited, timeDimID))
+  timeunit = 'seconds since '//adate(1:4)//'-'//adate(5:6)// &
+     '-'//adate(7:8)//' '//atime(1:2)//':'//atime(3:4)
+
+  ! lon
+     call nf90_err(nf90_def_dim(ncid, 'lon', nnx, lonDimID))
+  ! lat
+     call nf90_err(nf90_def_dim(ncid, 'lat', nny, latDimID))
+
+  ! level
+  call nf90_err(nf90_def_dim(ncid, 'level', numzgrid, levDimID))
+  ! number of species
+  call nf90_err(nf90_def_dim(ncid, 'numspec', nspec, nspecDimID))
+  ! number of release points
+  call nf90_err(nf90_def_dim(ncid, 'numpoint', numpoint, npointDimID))
+  ! number of age classes
+  call nf90_err(nf90_def_dim(ncid, 'nageclass', nageclass, nageclassDimID))
+  ! dimension for release point characters
+  call nf90_err(nf90_def_dim(ncid, 'nchar', 45, ncharDimID))
+
+  ! create variables
+  !*************************
+
+  ! time
+  call nf90_err(nf90_def_var(ncid, 'time', nf90_int, (/ timeDimID /), tID))
+  call nf90_err(nf90_put_att(ncid, tID, 'units', timeunit))
+  call nf90_err(nf90_put_att(ncid, tID, 'calendar', 'proleptic_gregorian'))
+  if (lnest) then
+     timeIDn = tID
+  else
+     timeID = tID
+  endif
+
+  ! lon
+    call nf90_err(nf90_def_var(ncid, 'lon', nf90_float, (/ lonDimID /), lonID))
+    call nf90_err(nf90_put_att(ncid, lonID, 'long_name', 'longitude in degree east'))
+
+  call nf90_err(nf90_put_att(ncid, lonID, 'axis', 'X'))
+  call nf90_err(nf90_put_att(ncid, lonID, 'units', 'degrees'))
+  call nf90_err(nf90_put_att(ncid, lonID, 'standard_name', 'grid_longitude'))
+  call nf90_err(nf90_put_att(ncid, lonID, 'description', 'grid cell centers'))
+  ! lat
+    call nf90_err(nf90_def_var(ncid, 'lat', nf90_float, (/ latDimID /), latID))
+    call nf90_err(nf90_put_att(ncid, latID, 'long_name', 'latitude in degree north'))
+
+  call nf90_err(nf90_put_att(ncid, latID, 'axis', 'Y'))
+  call nf90_err(nf90_put_att(ncid, latID, 'units', 'degrees'))
+  call nf90_err(nf90_put_att(ncid, latID, 'standard_name', 'grid_latitude'))
+  call nf90_err(nf90_put_att(ncid, latID, 'description', 'grid cell centers'))
+
+
+  ! level
+  call nf90_err(nf90_def_var(ncid, 'level', nf90_float, (/ levDimID /), levID))
+  call nf90_err(nf90_put_att(ncid, levID, 'axis', 'Z'))
+  call nf90_err(nf90_put_att(ncid, levID, 'units', 'meters'))
+  call nf90_err(nf90_put_att(ncid, levID, 'positive', 'up'))
+  call nf90_err(nf90_put_att(ncid, levID, 'standard_name', 'height'))
+  call nf90_err(nf90_put_att(ncid, levID, 'long_name', 'height above ground'))
+
+  ! release comment
+  call nf90_err(nf90_def_var(ncid, 'RELCOM', nf90_char, (/ ncharDimID,npointDimID /), &
+       relcomID))
+  call nf90_err(nf90_put_att(ncid, relcomID, 'long_name', 'release point name'))
+  ! release longitude 1
+  call nf90_err(nf90_def_var(ncid, 'RELLNG1', nf90_float, (/ npointDimID /), rellng1ID))
+  call nf90_err(nf90_put_att(ncid, rellng1ID, 'units', 'degrees_east'))
+  call nf90_err(nf90_put_att(ncid, rellng1ID, 'long_name', &
+       'release longitude lower left corner'))
+  ! release longitude 2
+  call nf90_err(nf90_def_var(ncid, 'RELLNG2', nf90_float, (/ npointDimID /), rellng2ID))
+  call nf90_err(nf90_put_att(ncid, rellng2ID, 'units', 'degrees_east'))
+  call nf90_err(nf90_put_att(ncid, rellng2ID, 'long_name', &
+       'release longitude upper right corner'))
+  ! release latitude 1
+  call nf90_err(nf90_def_var(ncid, 'RELLAT1', nf90_float, (/ npointDimID /), rellat1ID))
+  call nf90_err(nf90_put_att(ncid, rellat1ID, 'units', 'degrees_north'))
+  call nf90_err(nf90_put_att(ncid, rellat1ID, 'long_name', &
+       'release latitude lower left corner'))
+  ! release latitude 2
+  call nf90_err(nf90_def_var(ncid, 'RELLAT2', nf90_float, (/ npointDimID /), rellat2ID))
+  call nf90_err(nf90_put_att(ncid, rellat2ID, 'units', 'degrees_north'))
+  call nf90_err(nf90_put_att(ncid, rellat2ID, 'long_name', &
+       'release latitude upper right corner'))
+
+  ! hes: if rotated_ll it would be convenient also to write the the release points in rotated_coordinates
+
+  ! release height bottom
+  call nf90_err(nf90_def_var(ncid, 'RELZZ1', nf90_float, (/ npointDimID /), relzz1ID))
+  call nf90_err(nf90_put_att(ncid, relzz1ID, 'units', 'meters'))
+  call nf90_err(nf90_put_att(ncid, relzz1ID, 'long_name', 'release height bottom'))
+  ! release height top
+  call nf90_err(nf90_def_var(ncid, 'RELZZ2', nf90_float, (/ npointDimID /), relzz2ID))
+  call nf90_err(nf90_put_att(ncid, relzz2ID, 'units', 'meters'))
+  call nf90_err(nf90_put_att(ncid, relzz2ID, 'long_name', 'release height top'))
+  ! release kind
+  call nf90_err(nf90_def_var(ncid, 'RELKINDZ', nf90_int, (/ npointDimID /), relkindzID))
+  call nf90_err(nf90_put_att(ncid, relkindzID, 'long_name', 'release kind'))
+  ! release start
+  call nf90_err(nf90_def_var(ncid, 'RELSTART', nf90_int, (/ npointDimID /), relstartID))
+  call nf90_err(nf90_put_att(ncid, relstartID, 'units', 'seconds'))
+  call nf90_err(nf90_put_att(ncid, relstartID, 'long_name', &
+       'release start relative to simulation start'))
+  ! release end
+  call nf90_err(nf90_def_var(ncid, 'RELEND', nf90_int, (/ npointDimID /), relendID))
+  call nf90_err(nf90_put_att(ncid, relendID, 'units', 'seconds'))
+  call nf90_err(nf90_put_att(ncid, relendID, 'long_name', &
+       'release end relative to simulation start'))
+  ! release particles
+  call nf90_err(nf90_def_var(ncid, 'RELPART', nf90_int, (/ npointDimID /), relpartID))
+  call nf90_err(nf90_put_att(ncid, relpartID, 'long_name', 'number of release particles'))
+  ! release particle masses
+  call nf90_err(nf90_def_var(ncid, 'RELXMASS', nf90_float, (/ npointDimID, nspecDimID /), &
+       relxmassID))
+  call nf90_err(nf90_put_att(ncid, relxmassID, 'long_name', 'total release particle mass'))
+ 
+  ! age classes
+  call nf90_err(nf90_def_var(ncid, 'LAGE', nf90_int, (/ nageclassDimID /), lageID))
+  call nf90_err(nf90_put_att(ncid, lageID, 'units', 'seconds'))
+  call nf90_err(nf90_put_att(ncid, lageID, 'long_name', 'age class'))
+
+! output orography
+  if (.not. min_size) then
+    call nf90_err(nf90_def_var(ncid, 'ORO', nf90_int, (/ lonDimID, latDimID /), oroID,  &
+      deflate_level=deflate_level, chunksizes= (/ nnx, nny /)))
+    call nf90_err(nf90_put_att(ncid, oroID, 'standard_name', 'surface altitude'))
+    call nf90_err(nf90_put_att(ncid, oroID, 'long_name', 'outgrid surface altitude'))
+    call nf90_err(nf90_put_att(ncid, oroID, 'units', 'm'))
+  end if
+
+
+
+  ! concentration output, wet and dry deposition variables (one per species)
+  call output_units(units)
+
+  dIDs = (/ londimid, latdimid, levdimid, timedimid, npointdimid, nageclassdimid /)
+  depdIDs = (/ londimid, latdimid, timedimid, npointdimid, nageclassdimid /)
+  if (lnest) then
+     dimidsn    = dIDs
+     depdimidsn = depdIDs
+  else
+     dimids    = dIDs
+     depdimids = depdIDs
+  endif
+
+  ! set chunksizes according to largest written portion of data in an individual call to 
+  ! nf90_put_var
+  chunksizes = (/ nnx, nny, numzgrid, 1, 1, 1 /)
+  dep_chunksizes = (/ nnx, nny, 1, 1, 1 /)
+
+  do i = 1,nspec
+     write(anspec,'(i3.3)') i
+
+     ! concentration output
+     if (iout.eq.1.or.iout.eq.3.or.iout.eq.5) then
+        call nf90_err(nf90_def_var(ncid,'spec'//anspec, nf90_float, dIDs, sID , &
+             deflate_level = deflate_level,  &
+             chunksizes = chunksizes ))
+        call nf90_err(nf90_put_att(ncid, sID, 'units', units))
+        call nf90_err(nf90_put_att(ncid, sID, 'long_name', species(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'decay', decay(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'weightmolar', weightmolar(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
+        !call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)) ! what is this??
+        call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
+
+        if (lnest) then
+           specIDn(i) = sID
+        else
+           specID(i) = sID
+        endif
+     endif
+
+     ! mixing ratio output
+     if (iout.eq.2.or.iout.eq.3) then
+        call nf90_err(nf90_def_var(ncid,'spec'//anspec//'_pptv', nf90_float, dIDs, sID , &
+             deflate_level = deflate_level,  &
+             chunksizes = chunksizes ))
+        call nf90_err(nf90_put_att(ncid, sID, 'units', 'pptv'))
+        call nf90_err(nf90_put_att(ncid, sID, 'long_name', species(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'decay', decay(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'weightmolar', weightmolar(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'ohreact', ohreact(i)))
+        !call nf90_err(nf90_put_att(ncid, sID, 'kao', kao(i)) ! what is this??
+        call nf90_err(nf90_put_att(ncid, sID, 'vsetaver', vsetaver(i)))
+        call nf90_err(nf90_put_att(ncid, sID, 'spec_ass', spec_ass(i)))
+
+        if (lnest) then
+           specIDnppt(i) = sID
+        else
+           specIDppt(i) = sID
+        endif
+     endif
+
+     ! wet and dry deposition fields for forward runs
+     if (wetdep) then
+        call nf90_err(nf90_def_var(ncid,'WD_spec'//anspec, nf90_float, depdIDs, &
+             wdsID, deflate_level = deflate_level, &
+             chunksizes = dep_chunksizes))
+        call nf90_err(nf90_put_att(ncid, wdsID, 'units', '1e-12 kg m-2'))
+        call nf90_err(nf90_put_att(ncid, wdsID, 'weta', weta(i)))
+        call nf90_err(nf90_put_att(ncid, wdsID, 'wetb', wetb(i)))
+        call nf90_err(nf90_put_att(ncid, wdsID, 'dquer', dquer(i)))
+        call nf90_err(nf90_put_att(ncid, wdsID, 'henry', henry(i)))
+        if (lnest) then
+           wdspecIDn(i) = wdsID
+        else
+           wdspecID(i) = wdsID
+        endif
+     endif
+     if (drydep) then
+        call nf90_err(nf90_def_var(ncid,'DD_spec'//anspec, nf90_float, depdIDs, &
+             ddsID, deflate_level = deflate_level, &
+             chunksizes = dep_chunksizes))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'units', '1e-12 kg m-2'))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'dryvel', dryvel(i)))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'reldiff', reldiff(i)))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'henry', henry(i)))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'f0', f0(i)))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'dquer', dquer(i)))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'density', density(i)))
+        call nf90_err(nf90_put_att(ncid, ddsID, 'dsigma', dsigma(i)))
+        if (lnest) then
+           ddspecIDn(i) = ddsID
+        else
+           ddspecID(i) = ddsID
+        endif
+     endif
+  end do
+
+
+  ! global (metadata) attributes
+  !*******************************
+  call writemetadata(ncid,lnest)
+
+
+  ! moves the file from define to data mode
+  call nf90_err(nf90_enddef(ncid))
+
+!  ! hes: inquire var definition
+!  do i = 1,nspec
+!     write(anspec,'(i3.3)') i
+!
+!     ! concentration output
+!     if (iout.eq.1.or.iout.eq.3.or.iout.eq.5) then
+!        if (lnest) then
+!           sID = specIDn(i)
+!        else
+!           sID = specID(i)
+!        endif
+!        call nf90_err(nf90_inquire_variable(ncid, sID, chunksizes=inq_chunksizes))
+!        write(*,*) "Chunksizes for var "//anspec//": ", inq_chunksizes
+!     endif
+!  end do
+
+  
+  ! fill with data
+  !******************************
+  ! longitudes (grid cell centers)
+  if (lnest) then
+     do i = 1,numxgridn
+        coord(i) = outlon0n + (i-0.5)*dxoutn
+     enddo
+     call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgridn)))
+  else
+     do i = 1,numxgrid
+        coord(i) = outlon0 + (i-0.5)*dxout
+     enddo
+     call nf90_err(nf90_put_var(ncid, lonID, coord(1:numxgrid)))
+  endif
+  ! latitudes (grid cell centers)
+  if (lnest) then
+     do i = 1,numygridn
+        coord(i) = outlat0n + (i-0.5)*dyoutn
+     enddo
+     call nf90_err(nf90_put_var(ncid, latID, coord(1:numygridn)))
+  else
+     do i = 1,numygrid
+        coord(i) = outlat0 + (i-0.5)*dyout
+     enddo
+     call nf90_err(nf90_put_var(ncid, latID, coord(1:numygrid)))
+  endif
+  ! levels
+  call nf90_err(nf90_put_var(ncid, levID, outheight(1:numzgrid)))
+  
+  ! release point information
+  do i = 1,numpoint
+     call nf90_err(nf90_put_var(ncid, relstartID, ireleasestart(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, relendID, ireleaseend(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, relkindzID, kindz(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, rellng1ID, xpoint1(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, rellng2ID, xpoint2(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, rellat1ID, ypoint1(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, rellat2ID, ypoint2(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, relzz1ID, zpoint1(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, relzz2ID, zpoint2(i), (/i/)))
+     call nf90_err(nf90_put_var(ncid, relpartID, npart(i), (/i/)))
+     if (i .le. 1000) then
+       call nf90_err(nf90_put_var(ncid, relcomID, compoint(i), (/1,i/), (/45,1/)))
+     else
+       call nf90_err(nf90_put_var(ncid, relcomID, 'NA', (/1,i/), (/45,1/)))
+     endif 
+     call nf90_err(nf90_put_var(ncid, relxmassID, xmass(i,1:nspec), (/i,1/), (/1,nspec/)))
+  end do
+
+  ! age classes
+  call nf90_err(nf90_put_var(ncid, lageID, lage(1:nageclass)))
+
+  ! orography 
+  if (.not. min_size) then
+    if (lnest) then
+      call nf90_err(nf90_put_var(ncid, oroID, orooutn(0:(nnx-1), 0:(nny-1))))
+    else
+      call nf90_err(nf90_put_var(ncid, oroID, oroout(0:(nnx-1), 0:(nny-1))))
+    endif
+  end if
+
+  call nf90_err(nf90_close(ncid))
+
+  return
+
+end subroutine writeheader_ncdf
+
+
+subroutine concoutput_ncdf(itime,outnum,gridtotalunc,wetgridtotalunc, &
+     drygridtotalunc)
+  !                          i     i          o             o
+  !       o
+  !*****************************************************************************
+  !                                                                            *
+  !     Output of the concentration grid and the receptor concentrations.      *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     24 May 1995                                                            *
+  !                                                                            *
+  !     13 April 1999, Major update: if output size is smaller, dump output in *
+  !                    sparse matrix format; additional output of uncertainty  *
+  !                                                                            *
+  !     05 April 2000, Major update: output of age classes; output for backward*
+  !                    runs is time spent in grid cell times total mass of     *
+  !                    species.                                                *
+  !                                                                            *
+  !     17 February 2002, Appropriate dimensions for backward and forward runs *
+  !                       are now specified in module par_mod                  *
+  !                                                                            *
+  !     June 2006, write grid in sparse matrix with a single write command     *
+  !                in order to save disk space                                 *
+  !                                                                            *
+  !     2008 new sparse matrix format                                          *
+  !                                                                            *
+  !     February 2010, Dominik Brunner, Empa                                   *
+  !                    Adapted for COSMO                                       *
+  !                    Remark: calculation of density could be improved.       *
+  !                    Currently, it is calculated for the lower left corner   *
+  !                    of each output grid cell rather than for its center.    *
+  !                    Furthermore, the average density could be calculated    *
+  !                    from the difference in pressure at the top and bottom   *
+  !                    of each cell rather than by interpolation.              *
+  !                                                                            *
+  !     April 2013, Dominik Brunner, Empa                                      *
+  !                    Adapted for netcdf output                               *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! outnum          number of samples                                          *
+  ! ncells          number of cells with non-zero concentrations               *
+  ! sparse          .true. if in sparse matrix format, else .false.            *
+  ! tot_mu          1 for forward, initial mass mixing ration for backw. runs  *
+  !                                                                            *
+  !*****************************************************************************
+
+  use unc_mod, only: gridunc,drygridunc,wetgridunc
+
+
+  implicit none
+
+  integer, intent(in) :: itime
+  real, intent(in)    :: outnum
+  real, intent(out)   :: gridtotalunc,wetgridtotalunc,drygridtotalunc
+  real                :: densityoutrecept(maxreceptor)
+  integer             :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
+  integer             :: nage,i,l, jj
+  real                :: tot_mu(maxspec,maxpointspec_act)
+  real                :: halfheight,dz,dz1,dz2
+  real                :: xl,yl,xlrot,ylrot,zagnd,zagndprev
+  real                :: auxgrid(nclassunc),gridtotal,gridsigmatotal
+  real                :: wetgridtotal,wetgridsigmatotal
+  real                :: drygridtotal,drygridsigmatotal
+
+  real, parameter     :: weightair=28.97
+
+
+  ! open output file
+  call nf90_err(nf90_open(trim(ncfname), nf90_write, ncid))
+
+  ! write time
+  tpointer = tpointer + 1
+  call nf90_err(nf90_put_var( ncid, timeID, itime, (/ tpointer /)))
+  
+  ! For forward simulations, output fields have dimension MAXSPEC,
+  ! for backward simulations, output fields have dimension MAXPOINT.
+  ! Thus, make loops either about nspec, or about numpoint
+  !*****************************************************************
+
+  if (ldirect.eq.1) then
+    do ks=1,nspec
+      do kp=1,maxpointspec_act
+        tot_mu(ks,kp)=1.0
+      end do
+    end do
+  else
+    do ks=1,nspec
+      do kp=1,maxpointspec_act
+        tot_mu(ks,kp)=xmass(kp,ks)
+      end do
+    end do
+  endif
+
+
+  !*******************************************************************
+  ! Compute air density:
+  ! brd134: we now take into account whether we are in the mother or in
+  !    a nested domain (before only from mother domain)
+  ! Determine center altitude of output layer, and interpolate density
+  ! data to that altitude
+  !*******************************************************************
+
+  do kz=1,numzgrid
+    if (kz.eq.1) then
+      halfheight=outheight(1)/2.
+    else
+      halfheight=(outheight(kz)+outheight(kz-1))/2.
+    endif
+    do kzz=2,nz
+      if ((height(kzz-1).lt.halfheight).and. &
+           (height(kzz).gt.halfheight)) exit
+    end do
+    kzz=max(min(kzz,nz),2)
+    dz1=halfheight-height(kzz-1)
+    dz2=height(kzz)-halfheight
+    dz=dz1+dz2
+
+    do jy=0,numygrid-1
+      do ix=0,numxgrid-1
+        xl=outlon0+real(ix)*dxout
+        yl=outlat0+real(jy)*dyout
+        ! grid index in mother domain
+        xl=(xl-xlon0)/dx
+        yl=(yl-ylat0)/dx
+
+        ngrid=0
+        do jj=numbnests,1,-1
+          if ( xl.gt.xln(jj)+eps .and. xl.lt.xrn(jj)-eps .and. &
+                 yl.gt.yln(jj)+eps .and. yl.lt.yrn(jj)-eps ) then
+            ngrid=jj
+            exit 
+          end if
+        end do
+
+        if (ngrid.eq.0) then
+          iix=max(min(nint(xl),nxmin1),0) ! if output grid cell is outside mother domain
+          jjy=max(min(nint(yl),nymin1),0)
+
+          densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,memind(2))*dz1+ &
+             rho(iix,jjy,kzz-1,memind(2))*dz2)/dz
+        else
+          xl=(xl-xln(ngrid))*xresoln(ngrid)
+          yl=(yl-yln(ngrid))*yresoln(ngrid)
+          iix=max(min(nint(xl),nxn(ngrid)-1),0)
+          jjy=max(min(nint(yl),nyn(ngrid)-1),0)
+
+          densityoutgrid(ix,jy,kz)=(rhon(iix,jjy,kzz,memind(2), ngrid)*dz1+ &
+             rhon(iix,jjy,kzz-1,memind(2), ngrid)*dz2)/dz
+        endif
+      end do
+    end do
+  end do
+
+  ! brd134: for receptor points no option for nests yet to specify density
+  !    and also altitude zreceptor not considered yet (needs revision)
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    densityoutrecept(i)=rho(iix,jjy,1,memind(2))
+  end do
+
+  ! Output is different for forward and backward simulations
+  if (ldirect.eq.1) then
+     do kz=1,numzgrid
+        do jy=0,numygrid-1
+           do ix=0,numxgrid-1
+              factor3d(ix,jy,kz)=1.e12/volume(ix,jy,kz)/outnum
+           end do
+        end do
+     end do
+  else
+     do kz=1,numzgrid
+        do jy=0,numygrid-1
+           do ix=0,numxgrid-1
+              factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
+           end do
+        end do
+     end do
+  endif
+
+  !*********************************************************************
+  ! Determine the standard deviation of the mean concentration or mixing
+  ! ratio (uncertainty of the output) and the dry and wet deposition
+  !*********************************************************************
+
+  gridtotal=0.
+  gridsigmatotal=0.
+  gridtotalunc=0.
+  wetgridtotal=0.
+  wetgridsigmatotal=0.
+  wetgridtotalunc=0.
+  drygridtotal=0.
+  drygridsigmatotal=0.
+  drygridtotalunc=0.
+
+  do ks=1,nspec
+
+    do kp=1,maxpointspec_act
+      do nage=1,nageclass
+
+        do jy=0,numygrid-1
+          do ix=0,numxgrid-1
+
+            ! WET DEPOSITION
+            if ((wetdep).and.(ldirect.gt.0)) then
+              do l=1,nclassunc
+                auxgrid(l)=wetgridunc(ix,jy,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,wetgrid(ix,jy), &
+                   wetgridsigma(ix,jy),nclassunc)
+              ! Multiply by number of classes to get total concentration
+              wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
+              wetgridtotal=wetgridtotal+wetgrid(ix,jy)
+              ! Calculate standard deviation of the mean
+              wetgridsigma(ix,jy)= &
+                   wetgridsigma(ix,jy)* &
+                   sqrt(real(nclassunc))
+              wetgridsigmatotal=wetgridsigmatotal+ &
+                   wetgridsigma(ix,jy)
+            endif
+
+            ! DRY DEPOSITION
+            if ((drydep).and.(ldirect.gt.0)) then
+              do l=1,nclassunc
+                auxgrid(l)=drygridunc(ix,jy,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,drygrid(ix,jy), &
+                   drygridsigma(ix,jy),nclassunc)
+              ! Multiply by number of classes to get total concentration
+              drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc)
+              drygridtotal=drygridtotal+drygrid(ix,jy)
+              ! Calculate standard deviation of the mean
+              drygridsigma(ix,jy)= &
+                   drygridsigma(ix,jy)* &
+                   sqrt(real(nclassunc))
+              drygridsigmatotal=drygridsigmatotal+ &
+                   drygridsigma(ix,jy)
+            endif
+
+            ! CONCENTRATION OR MIXING RATIO
+            do kz=1,numzgrid
+              do l=1,nclassunc
+                auxgrid(l)=gridunc(ix,jy,kz,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,grid(ix,jy,kz), &
+                   gridsigma(ix,jy,kz),nclassunc)
+              ! Multiply by number of classes to get total concentration
+              grid(ix,jy,kz)= &
+                   grid(ix,jy,kz)*real(nclassunc)
+              gridtotal=gridtotal+grid(ix,jy,kz)
+              ! Calculate standard deviation of the mean
+              gridsigma(ix,jy,kz)= &
+                   gridsigma(ix,jy,kz)* &
+                   sqrt(real(nclassunc))
+              gridsigmatotal=gridsigmatotal+ &
+                   gridsigma(ix,jy,kz)
+            end do
+          end do
+        end do
+
+
+        !*******************************************************************
+        ! Generate output: may be in concentration (ng/m3) or in mixing
+        ! ratio (ppt) or both
+        ! Output the position and the values alternated multiplied by
+        ! 1 or -1, first line is number of values, number of positions
+        ! For backward simulations, the unit is seconds, stored in grid_time
+        !*******************************************************************
+
+        ! Concentration output
+        !*********************
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+
+          ! Wet deposition
+          if ((ldirect.eq.1).and.(WETDEP)) then
+             call nf90_err(nf90_put_var(ncid,wdspecID(ks),1.e12*&
+                  wetgrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
+          endif
+
+          ! Dry deposition
+          if ((ldirect.eq.1).and.(DRYDEP)) then
+             call nf90_err(nf90_put_var(ncid,ddspecID(ks),1.e12*&
+                  drygrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
+          endif
+
+          ! Concentrations
+          call nf90_err(nf90_put_var(ncid,specID(ks),grid(0:numxgrid-1,0:numygrid-1,&
+             1:numzgrid)*factor3d(0:numxgrid-1,0:numygrid-1,1:numzgrid)/tot_mu(ks,kp),&
+               (/ 1,1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /) ))
+ 
+        endif !  concentration output
+
+        ! Mixing ratio output
+        !********************
+
+        if ((iout.eq.2)) then      ! mixing ratio
+
+          ! Wet deposition
+          if ((ldirect.eq.1).and.(WETDEP)) then
+             call nf90_err(nf90_put_var(ncid,wdspecID(ks),1.e12*&
+                  wetgrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
+          endif
+
+          ! Dry deposition
+          if ((ldirect.eq.1).and.(DRYDEP)) then
+             call nf90_err(nf90_put_var(ncid,ddspecID(ks),1.e12*&
+                  drygrid(0:numxgrid-1,0:numygrid-1)/area(0:numxgrid-1,0:numygrid-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,1,1,1 /)))
+          endif
+
+          ! Mixing ratios
+          call nf90_err(nf90_put_var(ncid,specIDppt(ks),weightair/weightmolar(ks)*&
+               grid(0:numxgrid-1,0:numygrid-1,1:numzgrid)*&
+               factor3d(0:numxgrid-1,0:numygrid-1,1:numzgrid)/&
+               densityoutgrid(0:numxgrid-1,0:numygrid-1,1:numzgrid),&
+               (/ 1,1,1,tpointer,kp,nage /), (/ numxgrid,numygrid,numzgrid,1,1,1 /)))
+
+        endif ! output for ppt
+
+      end do
+    end do
+
+  end do
+
+  ! Close netCDF file
+  !**************************
+  call nf90_err(nf90_close(ncid))
+
+
+  if (gridtotal.gt.0.) gridtotalunc=gridsigmatotal/gridtotal
+  if (wetgridtotal.gt.0.) wetgridtotalunc=wetgridsigmatotal/ &
+       wetgridtotal
+  if (drygridtotal.gt.0.) drygridtotalunc=drygridsigmatotal/ &
+       drygridtotal
+
+  ! Dump of receptor concentrations
+
+  if (numreceptor.gt.0 .and. (iout.eq.2 .or. iout.eq.3)  ) then
+    write(unitoutreceptppt) itime
+    do ks=1,nspec
+      write(unitoutreceptppt) (1.e12*creceptor(i,ks)/outnum* &
+           weightair/weightmolar(ks)/densityoutrecept(i),i=1,numreceptor)
+    end do
+  endif
+
+  ! Dump of receptor concentrations
+
+  if (numreceptor.gt.0) then
+    write(unitoutrecept) itime
+    do ks=1,nspec
+      write(unitoutrecept) (1.e12*creceptor(i,ks)/outnum,i=1,numreceptor)
+    end do
+  endif
+
+
+  ! Reinitialization of grid
+  !*************************
+
+  creceptor(1:numreceptor,1:nspec) = 0.
+  gridunc(:,:,:,1:nspec,:,:,1:nageclass) = 0.  
+
+
+end subroutine concoutput_ncdf
+
+
+subroutine concoutput_nest_ncdf(itime,outnum)
+  !                               i     i 
+  !*****************************************************************************
+  !                                                                            *
+  !     Output of the concentration grid and the receptor concentrations.      *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     24 May 1995                                                            *
+  !                                                                            *
+  !     13 April 1999, Major update: if output size is smaller, dump output in *
+  !                    sparse matrix format; additional output of uncertainty  *
+  !                                                                            *
+  !     05 April 2000, Major update: output of age classes; output for backward*
+  !                    runs is time spent in grid cell times total mass of     *
+  !                    species.                                                *
+  !                                                                            *
+  !     17 February 2002, Appropriate dimensions for backward and forward runs *
+  !                    are now specified in module par_mod                     *
+  !                                                                            *
+  !     June 2006, write grid in sparse matrix with a single write command     *
+  !                    in order to save disk space                             *
+  !                                                                            *
+  !     2008 new sparse matrix format                                          *
+  !                                                                            *
+  !     19 February 2010, Dominik Brunner, Empa: Adapted for COSMO             *
+  !                                                                            *
+  !     April 2013, Dominik Brunner, Empa                                      *
+  !                    Adapted for netcdf output                               *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! itime           current simulation time                                    *
+  ! outnum          number of samples                                          *
+  !                                                                            *
+  !*****************************************************************************
+
+  use unc_mod, only: griduncn,drygriduncn,wetgriduncn
+ 
+  implicit none
+
+  integer, intent(in) :: itime
+  real, intent(in)    :: outnum
+  real                :: densityoutrecept(maxreceptor)
+  integer             :: ncid,kp,ks,kz,ix,jy,iix,jjy,kzz,kzzm1,ngrid
+  integer             :: nage,i,l, jj
+  real                :: tot_mu(maxspec,maxpointspec_act)
+  real                :: halfheight,dz,dz1,dz2
+  real                :: xl,yl,xlrot,ylrot,zagnd,zagndprev
+  real                :: auxgrid(nclassunc)
+  real, parameter     :: weightair=28.97
+
+  ! open output file
+  call nf90_err(nf90_open(trim(ncfnamen), nf90_write, ncid))
+
+  ! write time (do not increase time counter here, done in main output domain)
+  call nf90_err(nf90_put_var( ncid, timeID, itime, (/ tpointer /)))
+  
+  ! For forward simulations, output fields have dimension MAXSPEC,
+  ! for backward simulations, output fields have dimension MAXPOINT.
+  ! Thus, make loops either about nspec, or about numpoint
+  !*****************************************************************
+
+  if (ldirect.eq.1) then
+    do ks=1,nspec
+      do kp=1,maxpointspec_act
+        tot_mu(ks,kp)=1.0
+      end do
+    end do
+  else
+    do ks=1,nspec
+      do kp=1,maxpointspec_act
+        tot_mu(ks,kp)=xmass(kp,ks)
+      end do
+    end do
+  endif
+
+
+  !*******************************************************************
+  ! Compute air density:
+  ! brd134: we now take into account whether we are in the mother or in
+  !    a nested domain (before only from mother domain)
+  ! Determine center altitude of output layer, and interpolate density
+  ! data to that altitude
+  !*******************************************************************
+
+  do kz=1,numzgrid
+    if (kz.eq.1) then
+      halfheight=outheight(1)/2.
+    else
+      halfheight=(outheight(kz)+outheight(kz-1))/2.
+    endif
+    do kzz=2,nz
+      if ((height(kzz-1).lt.halfheight).and. &
+           (height(kzz).gt.halfheight)) exit
+    end do
+    kzz=max(min(kzz,nz),2)
+    dz1=halfheight-height(kzz-1)
+    dz2=height(kzz)-halfheight
+    dz=dz1+dz2
+
+    do jy=0,numygridn-1
+      do ix=0,numxgridn-1
+        xl=outlon0n+real(ix)*dxoutn
+        yl=outlat0n+real(jy)*dyoutn
+        xl=(xl-xlon0)/dx
+        yl=(yl-ylat0)/dy
+
+        ngrid=0
+        do jj=numbnests,1,-1
+          if ( xl.gt.xln(jj)+eps .and. xl.lt.xrn(jj)-eps .and. &
+                 yl.gt.yln(jj)+eps .and. yl.lt.yrn(jj)-eps ) then
+            ngrid=jj
+            exit 
+          end if
+        end do
+
+        if (ngrid.eq.0) then
+          iix=max(min(nint(xl),nxmin1),0)
+          jjy=max(min(nint(yl),nymin1),0)
+
+          densityoutgrid(ix,jy,kz)=(rho(iix,jjy,kzz,memind(2))*dz1+ &
+             rho(iix,jjy,kzz-1,memind(2))*dz2)/dz
+        else
+          xl=(xl-xln(ngrid))*xresoln(ngrid)
+          yl=(yl-yln(ngrid))*yresoln(ngrid)
+          iix=max(min(nint(xl),nxn(ngrid)-1),0)
+          jjy=max(min(nint(yl),nyn(ngrid)-1),0)
+          densityoutgrid(ix,jy,kz)=(rhon(iix,jjy,kzz,memind(2), ngrid)*dz1+ &
+             rhon(iix,jjy,kzz-1,memind(2), ngrid)*dz2)/dz
+        endif
+
+      end do
+    end do
+  end do
+
+  do i=1,numreceptor
+    xl=xreceptor(i)
+    yl=yreceptor(i)
+    iix=max(min(nint(xl),nxmin1),0)
+    jjy=max(min(nint(yl),nymin1),0)
+    densityoutrecept(i)=rho(iix,jjy,1,memind(2))
+  end do
+
+  ! Output is different for forward and backward simulations
+  if (ldirect.eq.1) then
+     do kz=1,numzgrid
+        do jy=0,numygridn-1
+           do ix=0,numxgridn-1
+              factor3d(ix,jy,kz)=1.e12/volumen(ix,jy,kz)/outnum
+           end do
+        end do
+     end do
+  else
+     do kz=1,numzgrid
+        do jy=0,numygridn-1
+           do ix=0,numxgridn-1
+              factor3d(ix,jy,kz)=real(abs(loutaver))/outnum
+           end do
+        end do
+     end do
+  endif
+
+  !*********************************************************************
+  ! Determine the standard deviation of the mean concentration or mixing
+  ! ratio (uncertainty of the output) and the dry and wet deposition
+  !*********************************************************************
+
+  do ks=1,nspec
+
+    do kp=1,maxpointspec_act
+      do nage=1,nageclass
+
+        do jy=0,numygridn-1
+          do ix=0,numxgridn-1
+            ! WET DEPOSITION
+            if ((WETDEP).and.(ldirect.gt.0)) then
+              do l=1,nclassunc
+                auxgrid(l)=wetgriduncn(ix,jy,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,wetgrid(ix,jy), &
+                   wetgridsigma(ix,jy),nclassunc)
+              ! Multiply by number of classes to get total concentration
+              wetgrid(ix,jy)=wetgrid(ix,jy)*real(nclassunc)
+              ! Calculate standard deviation of the mean
+              wetgridsigma(ix,jy)= &
+                   wetgridsigma(ix,jy)* &
+                   sqrt(real(nclassunc))
+            endif
+
+            ! DRY DEPOSITION
+            if ((DRYDEP).and.(ldirect.gt.0)) then
+              do l=1,nclassunc
+                auxgrid(l)=drygriduncn(ix,jy,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,drygrid(ix,jy), &
+                   drygridsigma(ix,jy),nclassunc)
+              ! Multiply by number of classes to get total concentration
+              drygrid(ix,jy)=drygrid(ix,jy)*real(nclassunc)
+              ! Calculate standard deviation of the mean
+              drygridsigma(ix,jy)= &
+                   drygridsigma(ix,jy)* &
+                   sqrt(real(nclassunc))
+            endif
+
+            ! CONCENTRATION OR MIXING RATIO
+            do kz=1,numzgrid
+              do l=1,nclassunc
+                auxgrid(l)=griduncn(ix,jy,kz,ks,kp,l,nage)
+              end do
+              call mean(auxgrid,grid(ix,jy,kz), &
+                   gridsigma(ix,jy,kz),nclassunc)
+              ! Multiply by number of classes to get total concentration
+              grid(ix,jy,kz)= &
+                   grid(ix,jy,kz)*real(nclassunc)
+              ! Calculate standard deviation of the mean
+              gridsigma(ix,jy,kz)= &
+                   gridsigma(ix,jy,kz)* &
+                   sqrt(real(nclassunc))
+            end do
+          end do
+        end do
+
+
+        !*******************************************************************
+        ! Generate output: may be in concentration (ng/m3) or in mixing
+        ! ratio (ppt) or both
+        ! Output the position and the values alternated multiplied by
+        ! 1 or -1, first line is number of values, number of positions
+        ! For backward simulations, the unit is seconds, stored in grid_time
+        !*******************************************************************
+
+        ! Concentration output
+        !*********************
+        if ((iout.eq.1).or.(iout.eq.3).or.(iout.eq.5)) then
+
+          ! Wet deposition
+          if ((ldirect.eq.1).and.(WETDEP)) then
+             call nf90_err(nf90_put_var(ncid,wdspecIDn(ks),1.e12*&
+                  wetgrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
+          endif
+
+          ! Dry deposition
+          if ((ldirect.eq.1).and.(DRYDEP)) then
+             call nf90_err(nf90_put_var(ncid,ddspecIDn(ks),1.e12*&
+                  drygrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
+          endif
+
+          ! Concentrations
+          call nf90_err(nf90_put_var(ncid,specIDn(ks),grid(0:numxgridn-1,0:numygridn-1,&
+             1:numzgrid)*factor3d(0:numxgridn-1,0:numygridn-1,1:numzgrid)/tot_mu(ks,kp),&
+               (/ 1,1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,numzgrid,1,1,1 /)))
+ 
+        endif !  concentration output
+
+        ! Mixing ratio output
+        !********************
+
+        if ((iout.eq.2)) then      ! mixing ratio
+
+          ! Wet deposition
+          if ((ldirect.eq.1).and.(WETDEP)) then
+             call nf90_err(nf90_put_var(ncid,wdspecIDn(ks),1.e12*&
+                  wetgrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
+          endif
+
+          ! Dry deposition
+          if ((ldirect.eq.1).and.(DRYDEP)) then
+             call nf90_err(nf90_put_var(ncid,ddspecIDn(ks),1.e12*&
+                  drygrid(0:numxgridn-1,0:numygridn-1)/arean(0:numxgridn-1,0:numygridn-1),&
+                  (/ 1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,1,1,1 /)))
+          endif
+
+          ! Mixing ratios
+          call nf90_err(nf90_put_var(ncid,specIDnppt(ks),weightair/weightmolar(ks)*&
+               grid(0:numxgridn-1,0:numygridn-1,1:numzgrid)*&
+               factor3d(0:numxgridn-1,0:numygridn-1,1:numzgrid)/&
+               densityoutgrid(0:numxgridn-1,0:numygridn-1,1:numzgrid),&
+               (/ 1,1,1,tpointer,kp,nage /), (/ numxgridn,numygridn,numzgrid,1,1,1 /)))
+
+        endif ! output for ppt
+
+      end do
+    end do
+
+  end do
+
+  ! Close netCDF file
+  !**************************
+  call nf90_err(nf90_close(ncid))
+
+
+  ! Reinitialization of grid
+  !*************************
+
+  creceptor(1:numreceptor,1:nspec) = 0.
+  griduncn(:,:,:,1:nspec,:,:,1:nageclass) = 0.  
+
+
+end subroutine concoutput_nest_ncdf
+
+
+end module nc_output_mod
+
+
+#endif
Index: /branches/flexpart91_hasod/src_parallel/timemanager.f90
===================================================================
--- /branches/flexpart91_hasod/src_parallel/timemanager.f90	(revision 8)
+++ /branches/flexpart91_hasod/src_parallel/timemanager.f90	(revision 9)
@@ -97,5 +97,5 @@
        idummy_ran3_initialize, idummy_ran3_advance
 #ifdef NETCDF_OUTPUT
-  use ncdf_output_mod, only: concoutput_ncdf, concoutput_nest_ncdf
+  use nc_output_mod, only: concoutput_ncdf, concoutput_nest_ncdf
 #endif 
 
