Index: src_parallel/FLEXPART.f90
===================================================================
--- src_parallel/FLEXPART.f90	(revision 7c1fd443e6ad66b1a84c706126014b304c08ff6e)
+++ src_parallel/FLEXPART.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -50,5 +50,5 @@
 
 #ifdef NETCDF_OUTPUT
-  use nc_output_mod, only: writeheader_ncdf
+  use netcdf_output_mod, only: writeheader_netcdf
 #endif
 
@@ -217,7 +217,7 @@
   !******************************************************************
 
-  if (lncdfout) then 
+  if (lnetcdfout.eq.1) then 
 #ifdef NETCDF_OUTPUT
-    call writeheader_ncdf(lnest = .false.)
+    call writeheader_netcdf(lnest = .false.)
 #endif
   else 
@@ -226,7 +226,7 @@
 
   if (nested_output.eq.1) then
-    if (lncdfout) then
+    if (lnetcdfout.eq.1) then
 #ifdef NETCDF_OUTPUT
-      call writeheader_ncdf(lnest = .true.)
+      call writeheader_netcdf(lnest = .true.)
 #endif
     else
@@ -235,5 +235,5 @@
   endif
 
-  if (.NOT.lncdfout) then
+  if (lnetcdfout.ne.1) then
     open(unitdates,file=path(2)(1:length(2))//'dates')
   end if
Index: src_parallel/Makefile.netcdf
===================================================================
--- src_parallel/Makefile.netcdf	(revision 7c1fd443e6ad66b1a84c706126014b304c08ff6e)
+++ src_parallel/Makefile.netcdf	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -25,5 +25,5 @@
 xmass_mod.o             flux_mod.o \
 point_mod.o             outg_mod.o \
-random_mod.o 		nc_output_mod.o
+random_mod.o 		netcdf_output_mod.o
 
 OBJECTS = \
Index: src_parallel/com_mod.f90
===================================================================
--- src_parallel/com_mod.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/com_mod.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -108,6 +108,6 @@
   ! lagespectra  1 if age spectra calculation switched on, 2 if not
 
-  logical :: lncdfout
-  !  switch for netcdf grid output, set to true if IOUT .gt. 10 in COMMAND (file/namelist)
+  integer :: lnetcdfout
+  ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
 
   integer :: nageclass,lage(maxageclass)
@@ -120,5 +120,4 @@
 
   ! gdomainfill             .T., if domain-filling is global, .F. if not
-
 
 
@@ -138,5 +137,4 @@
   !integer*2 kindz(maxpoint)
   integer :: specnum(maxspec)
-  !real xmass(maxpoint,maxspec)
   real :: decay(maxspec)
   real :: weta(maxspec),wetb(maxspec)
@@ -242,5 +240,5 @@
 
   integer :: numbwf,wftime(maxwf),lwindinterv
-  character(len=255) :: wfname(maxwf),wfspec(maxwf)
+  character(len=255) :: wfname(maxwf)
 
   ! lwindinterv [s]         Interval between wind fields currently in memory
@@ -248,6 +246,4 @@
   ! wftime(maxwf) [s]       times relative to beginning time of wind fields
   ! wfname(maxwf)           file names of wind fields
-  ! wfspec(maxwf)           specifications of wind field file, e.g. if on hard
-  !                    disc or on tape
 
   integer :: memtime(2),memind(2)
@@ -412,9 +408,6 @@
 
   character(len=255) :: wfnamen(maxnests,maxwf)
-  character(len=18) :: wfspecn(maxnests,maxwf)
 
   ! wfnamen      nested wind field names
-  ! wfspecn      specifications of wind field file, e.g. if on hard
-  !         disc or on tape
 
 
Index: src_parallel/com_mod.mod
===================================================================
--- src_parallel/com_mod.mod	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/com_mod.mod	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -1,4 +1,4 @@
-GFORTRAN module version '6' created from com_mod.f90 on Tue Jul 23 11:18:10 2013
-MD5:b4969303109781e9e60d59c7f3797b4e -- If you edit this, you'll get what you deserve.
+GFORTRAN module version '6' created from com_mod.f90 on Wed Aug 14 15:51:31 2013
+MD5:9175d0c3444530a2599b03a75393a9cb -- If you edit this, you'll get what you deserve.
 
 (() () () () () () () () () () () () () () () () () () () () () () ()
@@ -1217,37 +1217,26 @@
 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '5000')) 0 () () () 0 0)
-247 'wfspec' 'com_mod' 'wfspec' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (CHARACTER 1 0 0 CHARACTER ((
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '255'))) 0 0 () (1 0 EXPLICIT (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '5000')) 0 () () () 0 0)
-248 'wfspecn' 'com_mod' 'wfspecn' 1 ((VARIABLE UNKNOWN-INTENT
-UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (CHARACTER 1 0 0
-CHARACTER ((CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '18'))) 0 0 () (2 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '5000')) 0 () () () 0 0)
-249 'wftime' 'com_mod' 'wftime' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+247 'wftime' 'com_mod' 'wftime' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (INTEGER 4 0 0 INTEGER ()) 0 0 () (
 1 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (
 INTEGER 4 0 0 INTEGER ()) 0 '5000')) 0 () () () 0 0)
-250 'wstar' 'com_mod' 'wstar' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (4 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '360') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '180') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '2')) 0 () () () 0 0)
-251 'wstarn' 'com_mod' 'wstarn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (5 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '2') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-252 'ww' 'com_mod' 'ww' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+248 'wstar' 'com_mod' 'wstar' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (4 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '360') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '180') (CONSTANT (INTEGER 4 0 0
+INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
+INTEGER ()) 0 '2')) 0 () () () 0 0)
+249 'wstarn' 'com_mod' 'wstarn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (5 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0
+INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
+INTEGER ()) 0 '2') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+250 'ww' 'com_mod' 'ww' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
 IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (4 0 EXPLICIT (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER 4 0 0
@@ -1257,17 +1246,17 @@
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
 INTEGER ()) 0 '2')) 0 () () () 0 0)
-253 'wwn' 'com_mod' 'wwn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (5 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '92') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
-INTEGER ()) 0 '2') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
-CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-254 'xglobal' 'com_mod' 'xglobal' 1 ((VARIABLE UNKNOWN-INTENT
+251 'wwn' 'com_mod' 'wwn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (5 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '-1') (CONSTANT (INTEGER 4 0 0
+INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '92') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
+INTEGER ()) 0 '2') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
+CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+252 'xglobal' 'com_mod' 'xglobal' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0) (LOGICAL 4 0 0 LOGICAL ()) 0 0 ()
 () 0 () () () 0 0)
-255 'xlanduse' 'com_mod' 'xlanduse' 1 ((VARIABLE UNKNOWN-INTENT
+253 'xlanduse' 'com_mod' 'xlanduse' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (3 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
@@ -1276,5 +1265,5 @@
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
 INTEGER ()) 0 '13')) 0 () () () 0 0)
-256 'xlandusen' 'com_mod' 'xlandusen' 1 ((VARIABLE UNKNOWN-INTENT
+254 'xlandusen' 'com_mod' 'xlandusen' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (4 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
@@ -1284,81 +1273,81 @@
 INTEGER ()) 0 '13') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-257 'xln' 'com_mod' 'xln' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-258 'xlon0' 'com_mod' 'xlon0' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+255 'xln' 'com_mod' 'xln' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+256 'xlon0' 'com_mod' 'xlon0' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-259 'xlon0n' 'com_mod' 'xlon0n' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-260 'xmass1' 'com_mod' 'xmass1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+257 'xlon0n' 'com_mod' 'xlon0n' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+258 'xmass1' 'com_mod' 'xmass1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (2 0
 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
 4 0 0 INTEGER ()) 0 '3500000') (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')
 (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-261 'xmassperparticle' 'com_mod' 'xmassperparticle' 1 ((VARIABLE
+259 'xmassperparticle' 'com_mod' 'xmassperparticle' 1 ((VARIABLE
 UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL
 ()) 0 0 () () 0 () () () 0 0)
-262 'xoutshift' 'com_mod' 'xoutshift' 1 ((VARIABLE UNKNOWN-INTENT
+260 'xoutshift' 'com_mod' 'xoutshift' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL ()) 0 0 () () 0
 () () () 0 0)
-263 'xoutshiftn' 'com_mod' 'xoutshiftn' 1 ((VARIABLE UNKNOWN-INTENT
+261 'xoutshiftn' 'com_mod' 'xoutshiftn' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL ()) 0 0 () () 0
 () () () 0 0)
-264 'xreceptor' 'com_mod' 'xreceptor' 1 ((VARIABLE UNKNOWN-INTENT
+262 'xreceptor' 'com_mod' 'xreceptor' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (1 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '20')) 0 () () () 0 0)
-265 'xresoln' 'com_mod' 'xresoln' 1 ((VARIABLE UNKNOWN-INTENT
+263 'xresoln' 'com_mod' 'xresoln' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (1 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-266 'xrn' 'com_mod' 'xrn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-267 'xtra1' 'com_mod' 'xtra1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+264 'xrn' 'com_mod' 'xrn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+265 'xtra1' 'com_mod' 'xtra1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 8 0 0 REAL ()) 0 0 () (1 0
 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
 4 0 0 INTEGER ()) 0 '3500000')) 0 () () () 0 0)
-268 'ylat0' 'com_mod' 'ylat0' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+266 'ylat0' 'com_mod' 'ylat0' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL ()) 0 0 () () 0 () () () 0 0)
-269 'ylat0n' 'com_mod' 'ylat0n' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-270 'yln' 'com_mod' 'yln' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-271 'youtshift' 'com_mod' 'youtshift' 1 ((VARIABLE UNKNOWN-INTENT
+267 'ylat0n' 'com_mod' 'ylat0n' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+268 'yln' 'com_mod' 'yln' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+269 'youtshift' 'com_mod' 'youtshift' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL ()) 0 0 () () 0
 () () () 0 0)
-272 'youtshiftn' 'com_mod' 'youtshiftn' 1 ((VARIABLE UNKNOWN-INTENT
+270 'youtshiftn' 'com_mod' 'youtshiftn' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0) (REAL 4 0 0 REAL ()) 0 0 () () 0
 () () () 0 0)
-273 'yreceptor' 'com_mod' 'yreceptor' 1 ((VARIABLE UNKNOWN-INTENT
+271 'yreceptor' 'com_mod' 'yreceptor' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (1 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '20')) 0 () () () 0 0)
-274 'yresoln' 'com_mod' 'yresoln' 1 ((VARIABLE UNKNOWN-INTENT
+272 'yresoln' 'com_mod' 'yresoln' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (1 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '0') (
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-275 'yrn' 'com_mod' 'yrn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
-UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
-EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
-4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
-276 'ytra1' 'com_mod' 'ytra1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+273 'yrn' 'com_mod' 'yrn' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
+EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
+4 0 0 INTEGER ()) 0 '1')) 0 () () () 0 0)
+274 'ytra1' 'com_mod' 'ytra1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 8 0 0 REAL ()) 0 0 () (1 0
 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
 4 0 0 INTEGER ()) 0 '3500000')) 0 () () () 0 0)
-277 'z0' 'com_mod' 'z0' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
+275 'z0' 'com_mod' 'z0' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC UNKNOWN
 IMPLICIT-SAVE 0 0 DIMENSION THREADPRIVATE) (REAL 4 0 0 REAL ()) 0 0 () (
 1 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (
 INTEGER 4 0 0 INTEGER ()) 0 '13')) 0 () () () 0 0)
-278 'zcolumn_sn' 'com_mod' 'zcolumn_sn' 1 ((VARIABLE UNKNOWN-INTENT
+276 'zcolumn_sn' 'com_mod' 'zcolumn_sn' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (3 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
@@ -1367,5 +1356,5 @@
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
 INTEGER ()) 0 '3000')) 0 () () () 0 0)
-279 'zcolumn_we' 'com_mod' 'zcolumn_we' 1 ((VARIABLE UNKNOWN-INTENT
+277 'zcolumn_we' 'com_mod' 'zcolumn_we' 1 ((VARIABLE UNKNOWN-INTENT
 UNKNOWN-PROC UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0
 0 () (3 0 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (
@@ -1374,5 +1363,5 @@
 CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER 4 0 0
 INTEGER ()) 0 '3000')) 0 () () () 0 0)
-280 'ztra1' 'com_mod' 'ztra1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
+278 'ztra1' 'com_mod' 'ztra1' 1 ((VARIABLE UNKNOWN-INTENT UNKNOWN-PROC
 UNKNOWN IMPLICIT-SAVE 0 0 DIMENSION) (REAL 4 0 0 REAL ()) 0 0 () (1 0
 EXPLICIT (CONSTANT (INTEGER 4 0 0 INTEGER ()) 0 '1') (CONSTANT (INTEGER
@@ -1430,10 +1419,10 @@
 0 234 'vdepn' 0 235 'vset' 0 236 'vsetaver' 0 237 'vv' 0 238 'vvn' 0 239
 'vvpol' 0 240 'weightmolar' 0 241 'weta' 0 242 'wetb' 0 243 'wetdep' 0
-244 'wfname' 0 245 'wfnamen' 0 246 'wfspec' 0 247 'wfspecn' 0 248 'wftime'
-0 249 'wstar' 0 250 'wstarn' 0 251 'ww' 0 252 'wwn' 0 253 'xglobal' 0
-254 'xlanduse' 0 255 'xlandusen' 0 256 'xln' 0 257 'xlon0' 0 258 'xlon0n'
-0 259 'xmass1' 0 260 'xmassperparticle' 0 261 'xoutshift' 0 262
-'xoutshiftn' 0 263 'xreceptor' 0 264 'xresoln' 0 265 'xrn' 0 266 'xtra1'
-0 267 'ylat0' 0 268 'ylat0n' 0 269 'yln' 0 270 'youtshift' 0 271
-'youtshiftn' 0 272 'yreceptor' 0 273 'yresoln' 0 274 'yrn' 0 275 'ytra1'
-0 276 'z0' 0 277 'zcolumn_sn' 0 278 'zcolumn_we' 0 279 'ztra1' 0 280)
+244 'wfname' 0 245 'wfnamen' 0 246 'wftime' 0 247 'wstar' 0 248 'wstarn'
+0 249 'ww' 0 250 'wwn' 0 251 'xglobal' 0 252 'xlanduse' 0 253 'xlandusen'
+0 254 'xln' 0 255 'xlon0' 0 256 'xlon0n' 0 257 'xmass1' 0 258
+'xmassperparticle' 0 259 'xoutshift' 0 260 'xoutshiftn' 0 261 'xreceptor'
+0 262 'xresoln' 0 263 'xrn' 0 264 'xtra1' 0 265 'ylat0' 0 266 'ylat0n' 0
+267 'yln' 0 268 'youtshift' 0 269 'youtshiftn' 0 270 'yreceptor' 0 271
+'yresoln' 0 272 'yrn' 0 273 'ytra1' 0 274 'z0' 0 275 'zcolumn_sn' 0 276
+'zcolumn_we' 0 277 'ztra1' 0 278)
Index: src_parallel/nc_output_mod.f90
===================================================================
--- src_parallel/nc_output_mod.f90	(revision 7c1fd443e6ad66b1a84c706126014b304c08ff6e)
+++ 	(revision )
@@ -1,1295 +1,0 @@
-!**********************************************************************
-! 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: src_parallel/netcdf_output_mod.f90
===================================================================
--- src_parallel/netcdf_output_mod.f90	(revision None)
+++ src_parallel/netcdf_output_mod.f90	(revision None)
@@ -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_netcdf generates file including all information previously    *
+  !    stored in separate header files                                         *
+  !  - concoutput_netcdf write concentration output and wet and dry deposition   *
+  !                                                                            *
+  !     Author: D. Brunner                                                     *
+  !                                                                            *
+  !     12 April 2013                                                          *
+  !                                                                            *
+  !*****************************************************************************
+
+#ifdef NETCDF_OUTPUT
+
+module netcdf_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_netcdf and concoutput_nest_netcdf
+                       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_netcdf(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_netcdf
+
+
+subroutine concoutput_netcdf(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_netcdf
+
+
+subroutine concoutput_nest_netcdf(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_netcdf
+
+
+end module netcdf_output_mod
+
+
+#endif
Index: src_parallel/readageclasses.f90
===================================================================
--- src_parallel/readageclasses.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readageclasses.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -28,6 +28,7 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
-  !                                                                            *
   !     20 March 2000                                                          *
+  !     HSO, 14 August 2013
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -46,9 +47,17 @@
   integer :: i
 
+  ! namelist help variables
+  integer :: readerror
+
+  ! namelist declaration
+  namelist /ageclass/ &
+    nageclass, &
+    lage
+
+  nageclass=-1 ! preset to negative value to identify failed namelist input
 
   ! If age spectra calculation is switched off, set number of age classes
   ! to 1 and maximum age to a large number
   !**********************************************************************
-
 
   if (lagespectra.ne.1) then
@@ -57,5 +66,4 @@
     return
   endif
-
 
   ! If age spectra claculation is switched on,
@@ -66,9 +74,20 @@
        status='old',err=999)
 
-  do i=1,13
-    read(unitageclasses,*)
-  end do
-  read(unitageclasses,*) nageclass
+  ! try to read in as a namelist
+  read(unitageclasses,ageclass,iostat=readerror)
 
+  if ((nageclass.lt.0).or.(readerror.ne.0)) then
+    rewind(unitageclasses)
+    do i=1,13
+      read(unitageclasses,*)
+    end do
+    read(unitageclasses,*) nageclass
+    read(unitageclasses,*) lage(1)
+    do i=2,nageclass
+      read(unitageclasses,*) lage(i)
+    end do
+  endif
+
+  close(unitageclasses)
 
   if (nageclass.gt.maxageclass) then
@@ -81,5 +100,4 @@
   endif
 
-  read(unitageclasses,*) lage(1)
   if (lage(1).le.0) then
     write(*,*) ' #### FLEXPART MODEL ERROR! AGE OF FIRST      #### '
@@ -90,5 +108,4 @@
 
   do i=2,nageclass
-    read(unitageclasses,*) lage(i)
     if (lage(i).le.lage(i-1)) then
       write(*,*) ' #### FLEXPART MODEL ERROR! AGE CLASSES     #### '
Index: src_parallel/readavailable.f90
===================================================================
--- src_parallel/readavailable.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readavailable.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -46,7 +46,6 @@
   ! numbwf               actual number of wind fields                          *
   ! wfname(maxwf)        file names of needed wind fields                      *
-  ! wfspec(maxwf)        file specifications of wind fields (e.g., if on disc) *
   ! wftime(maxwf) [s]times of wind fields relative to beginning time           *
-  ! wfname1,wfspec1,wftime1 = same as above, but only local (help variables)   *
+  ! wfname1,wftime1 = same as above, but only local (help variables)   *
   !                                                                            *
   ! Constants:                                                                 *
@@ -64,7 +63,6 @@
   integer :: wftime1n(maxnests,maxwf),wftimen(maxnests,maxwf)
   real(kind=dp) :: juldate,jul,beg,end
-  character(len=255) :: fname,spec,wfname1(maxwf),wfspec1(maxwf)
+  character(len=255) :: fname,spec,wfname1(maxwf)
   character(len=255) :: wfname1n(maxnests,maxwf)
-  character(len=40) :: wfspec1n(maxnests,maxwf)
 
 
@@ -109,5 +107,4 @@
 
       wfname1(numbwf)=fname(1:index(fname,' '))
-      wfspec1(numbwf)=spec
       wftime1(numbwf)=nint((jul-bdate)*86400._dp)
     endif
@@ -146,5 +143,4 @@
 
         wfname1n(k,numbwfn(k))=fname
-        wfspec1n(k,numbwfn(k))=spec
         wftime1n(k,numbwfn(k))=nint((jul-bdate)*86400._dp)
       endif
@@ -205,5 +201,4 @@
     do i=1,numbwf
       wfname(i)=wfname1(i)
-      wfspec(i)=wfspec1(i)
       wftime(i)=wftime1(i)
     end do
@@ -211,5 +206,4 @@
       do i=1,numbwfn(k)
         wfnamen(k,i)=wfname1n(k,i)
-        wfspecn(k,i)=wfspec1n(k,i)
         wftimen(k,i)=wftime1n(k,i)
       end do
@@ -218,5 +212,4 @@
     do i=1,numbwf
       wfname(numbwf-i+1)=wfname1(i)
-      wfspec(numbwf-i+1)=wfspec1(i)
       wftime(numbwf-i+1)=wftime1(i)
     end do
@@ -224,5 +217,4 @@
       do i=1,numbwfn(k)
         wfnamen(k,numbwfn(k)-i+1)=wfname1n(k,i)
-        wfspecn(k,numbwfn(k)-i+1)=wfspec1n(k,i)
         wftimen(k,numbwfn(k)-i+1)=wftime1n(k,i)
       end do
Index: src_parallel/readcommand.f90
===================================================================
--- src_parallel/readcommand.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readcommand.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -27,6 +27,7 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
-  !                                                                            *
   !     18 May 1996                                                            *
+  !     HSO, 14 August 2013
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -83,4 +84,5 @@
   integer :: readerror
 
+  ! declaration of namelist
   namelist /command/ &
     ldirect, &
@@ -96,4 +98,5 @@
     iout, &
     ipout, &
+    lnetcdfout, &
     lsubgrid, &
     lconvection, &
@@ -109,9 +112,9 @@
     linit_cond
 
-  ! Presetting namelist command
+  ! Presetting the namelist command
   ldirect=0
-  ibdate=20000101
+  ibdate=0
   ibtime=0
-  iedate=20000102
+  iedate=0
   ietime=0
   loutstep=10800
@@ -124,4 +127,5 @@
   iout=3
   ipout=0
+  lnetcdfout=0
   lsubgrid=1
   lconvection=1
@@ -141,22 +145,14 @@
   !**************************************************************************
   open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
-         form='formatted',iostat=readerror)
-  ! If fail, check if file does not exist
-  if (readerror.ne.0) then
-
-    print*,'***ERROR: file COMMAND not found. Check your pathnames file.'
-    stop
-
-  endif
-
+         err=999)
+
+  ! try namelist input
   read(unitcommand,command,iostat=readerror)
-  close(unitcommand)
 
   ! If error in namelist format, try to open with old input code
   ! Failsafe detection with ldirect initial value change
-  if ((readerror.ne.0).or.(ldirect.ne.0)) then
-
-    open(unitcommand,file=path(1)(1:length(1))//'COMMAND',status='old', &
-         err=999)
+  if ((readerror.ne.0).or.(ldirect.eq.0)) then
+
+    rewind(unitcommand)
 
     ! Check the format of the COMMAND file (either in free format,
@@ -346,11 +342,10 @@
   endif
 
-!  check for netcdf output switch
-  lncdfout = .false.
+!  check for netcdf output switch (use for non-namelist input only!)
   if (iout.ge.8) then
-     lncdfout = .true.
+     lnetcdfout = 1
      iout = iout - 8
 #ifndef NETCDF_OUTPUT
-     print*,'ERROR: netcdf output not activated during compile time'
+     print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
 #endif
   endif
Index: src_parallel/readoutgrid.f90
===================================================================
--- src_parallel/readoutgrid.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readoutgrid.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -53,5 +53,19 @@
   real,parameter :: eps=1.e-4
 
-
+  ! namelist variables
+  integer, parameter :: maxoutlev=500
+  real :: outheights(maxoutlev)
+  integer :: readerror
+
+  ! declare namelist
+  namelist /outgrid/ &
+    outlon0,outlat0, &
+    numxgrid,numygrid, &
+    dxout,dyout, &
+    outheights
+
+  ! helps identifying failed namelist input
+  dxout=-1.0
+  outheights=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
@@ -61,24 +75,32 @@
        err=999)
 
-
-  call skplin(5,unitoutgrid)
-
-
-  ! 1.  Read horizontal grid specifications
-  !****************************************
-
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlon0
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlat0
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numxgrid
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numygrid
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dxout
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dyout
-
+  ! try namelist input
+  read(unitoutgrid,outgrid,iostat=readerror)
+
+
+  if ((dxout.le.0).or.(readerror.ne.0)) then
+
+    readerror=1
+
+    rewind(unitoutgrid)
+    call skplin(5,unitoutgrid)
+
+    ! 1.  Read horizontal grid specifications
+    !****************************************
+
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlon0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlat0
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numxgrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numygrid
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dxout
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dyout
+
+  endif
 
   ! Check validity of output grid (shall be within model domain)
@@ -102,6 +124,8 @@
   ! 2. Count Vertical levels of output grid
   !****************************************
-  j=0
-100   j=j+1
+
+  if (readerror.ne.0) then
+    j=0
+100 j=j+1
     do i=1,3
       read(unitoutgrid,*,end=99)
@@ -110,31 +134,42 @@
     if (outhelp.eq.0.) goto 99
     goto 100
-99   numzgrid=j-1
-
-    allocate(outheight(numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(outheighthalf(numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-
-
-  rewind(unitoutgrid)
-  call skplin(29,unitoutgrid)
+99  numzgrid=j-1
+  else
+    do i=1,maxoutlev
+      if (outheights(i).lt.0) exit
+    end do
+    numzgrid=i-1
+  end if
+
+  allocate(outheight(numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheight'
+  allocate(outheighthalf(numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate outheighthalf'
 
   ! 2. Vertical levels of output grid
   !**********************************
 
-  j=0
-1000   j=j+1
-    do i=1,3
-      read(unitoutgrid,*,end=990)
-    end do
-    read(unitoutgrid,'(4x,f7.1)',end=990) outhelp
-    if (outhelp.eq.0.) goto 99
-    outheight(j)=outhelp
-    goto 1000
-990   numzgrid=j-1
-
+  if (readerror.ne.0) then
+
+    rewind(unitoutgrid)
+    call skplin(29,unitoutgrid)
+
+    do j=1,numzgrid
+      do i=1,3
+        read(unitoutgrid,*)
+      end do
+      read(unitoutgrid,'(4x,f7.1)') outhelp
+      outheight(j)=outhelp
+    end do
+
+  else
+
+    do j=1,numzgrid
+      outheight(j)=outheights(j)
+    end do
+
+  endif
+
+  close(unitoutgrid)
 
   ! Check whether vertical levels are specified in ascending order
@@ -158,28 +193,21 @@
   end do
 
-
   xoutshift=xlon0-outlon0
   youtshift=ylat0-outlat0
-  close(unitoutgrid)
-
-    allocate(oroout(0:numxgrid-1,0:numygrid-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(area(0:numxgrid-1,0:numygrid-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+
+  allocate(oroout(0:numxgrid-1,0:numygrid-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate oroout'
+  allocate(area(0:numxgrid-1,0:numygrid-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate area'
+  allocate(volume(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate volume'
+  allocate(areaeast(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate areaeast'
+  allocate(areanorth(0:numxgrid-1,0:numygrid-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate areanorth'
   return
 
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "OUTGRID"    #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
   write(*,*) ' #### xxx/flexpart/options                    #### '
Index: src_parallel/readoutgrid_nest.f90
===================================================================
--- src_parallel/readoutgrid_nest.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readoutgrid_nest.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -53,5 +53,14 @@
   real,parameter :: eps=1.e-4
 
+  integer :: readerror
 
+  ! declare namelist
+  namelist /outgridn/ &
+    outlon0n,outlat0n, &
+    numxgridn,numygridn, &
+    dxoutn,dyoutn
+
+  ! helps identifying failed namelist input
+  dxoutn=-1.0
 
   ! Open the OUTGRID file and read output grid specifications
@@ -61,34 +70,38 @@
        status='old',err=999)
 
+  ! try namelist input
+  read(unitoutgrid,outgridn,iostat=readerror)
 
-  call skplin(5,unitoutgrid)
+  if ((dxoutn.le.0).or.(readerror.ne.0)) then
 
+    rewind(unitoutgrid)
+    call skplin(5,unitoutgrid)
 
-  ! 1.  Read horizontal grid specifications
-  !****************************************
+    ! 1.  Read horizontal grid specifications
+    !****************************************
 
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlon0n
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f11.4)') outlat0n
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numxgridn
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,i5)') numygridn
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dxoutn
-  call skplin(3,unitoutgrid)
-  read(unitoutgrid,'(4x,f12.5)') dyoutn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlon0n
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f11.4)') outlat0n
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numxgridn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,i5)') numygridn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dxoutn
+    call skplin(3,unitoutgrid)
+    read(unitoutgrid,'(4x,f12.5)') dyoutn
 
+  endif
 
-    allocate(orooutn(0:numxgridn-1,0:numygridn-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(arean(0:numxgridn-1,0:numygridn-1) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
-    allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate outh'
+  close(unitoutgrid)
+
+  allocate(orooutn(0:numxgridn-1,0:numygridn-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate orooutn'
+  allocate(arean(0:numxgridn-1,0:numygridn-1),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate arean'
+  allocate(volumen(0:numxgridn-1,0:numygridn-1,numzgrid),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate volumen'
 
   ! Check validity of output grid (shall be within model domain)
@@ -110,9 +123,7 @@
   xoutshiftn=xlon0-outlon0n
   youtshiftn=ylat0-outlat0n
-  close(unitoutgrid)
   return
 
-
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE OUTGRID_NEST #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
   write(*,*) ' #### xxx/flexpart/options                    #### '
Index: src_parallel/readreceptors.f90
===================================================================
--- src_parallel/readreceptors.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readreceptors.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -27,6 +27,7 @@
   !                                                                            *
   !     Author: A. Stohl                                                       *
-  !                                                                            *
   !     1 August 1996                                                          *
+  !     HSO, 14 August 2013
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -51,4 +52,12 @@
   character(len=16) :: receptor
 
+  integer :: readerror
+  real :: lon,lat   ! for namelist input, lon/lat are used instead of x,y
+
+  ! declare namelist
+  namelist /receptors/ &
+    receptor, lon, lat
+
+  lon=-999.9
 
   ! For backward runs, do not allow receptor output. Thus, set number of receptors to zero
@@ -67,12 +76,17 @@
        status='old',err=999)
 
-  call skplin(5,unitreceptor)
+  ! try namelist input
+  read(unitreceptor,receptors,iostat=readerror)
 
+  rewind(unitreceptor)
+  if ((lon.lt.-900).or.(readerror.ne.0)) then
 
-  ! Read the names and coordinates of the receptors
-  !************************************************
+    call skplin(5,unitreceptor)
 
-  j=0
-100   j=j+1
+    ! Read the names and coordinates of the receptors
+    !************************************************
+
+    j=0
+100 j=j+1
     read(unitreceptor,*,end=99)
     read(unitreceptor,*,end=99)
@@ -89,8 +103,8 @@
     endif
     if (j.gt.maxreceptor) then
-    write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
-    write(*,*) ' #### POINTS ARE GIVEN.                       #### '
-    write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
-    write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
+      write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
+      write(*,*) ' #### POINTS ARE GIVEN.                       #### '
+      write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
+      write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
     endif
     receptorname(j)=receptor
@@ -102,10 +116,40 @@
     goto 100
 
-99   numreceptor=j-1
+99  numreceptor=j-1
+
+  else ! continue with namelist input
+
+    j=0
+    do while (readerror.eq.0)
+      j=j+1
+      lon=-999.9
+      read(unitreceptor,receptors,iostat=readerror)
+      if ((lon.lt.-900).or.(readerror.ne.0)) then
+        readerror=1
+      else
+        if (j.gt.maxreceptor) then
+          write(*,*) ' #### FLEXPART MODEL ERROR! TOO MANY RECEPTOR #### '
+          write(*,*) ' #### POINTS ARE GIVEN.                       #### '
+          write(*,*) ' #### MAXIMUM NUMBER IS ',maxreceptor,'       #### '
+          write(*,*) ' #### PLEASE MAKE CHANGES IN FILE RECEPTORS   #### '
+        endif
+        receptorname(j)=receptor
+        xreceptor(j)=(lon-xlon0)/dx       ! transform to grid coordinates
+        yreceptor(j)=(lat-ylat0)/dy
+        xm=r_earth*cos(lat*pi/180.)*dx/180.*pi
+        ym=r_earth*dy/180.*pi
+        receptorarea(j)=xm*ym
+      endif
+    end do
+    numreceptor=j-1
+
+  endif
+
+
   close(unitreceptor)
   return
 
 
-999   write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
+999 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"  #### '
   write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
   write(*,'(a)') path(1)(1:length(1))
Index: src_parallel/readreleases.f90
===================================================================
--- src_parallel/readreleases.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readreleases.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -33,4 +33,6 @@
   !     Update: 29 January 2001                                                *
   !     Release altitude can be either in magl or masl                         *
+  !     HSO, 12 August 2013
+  !     Added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -67,5 +69,6 @@
   implicit none
 
-  integer :: numpartmax,i,j,id1,it1,id2,it2,specnum_rel,idum,stat
+  integer :: numpartmax,i,j,id1,it1,id2,it2,idum,stat
+  integer :: specnum_rel(maxspec)
   real :: vsh(ni),fracth(ni),schmih(ni),releaserate,xdum,cun
   real(kind=dp) :: jul1,jul2,juldate
@@ -73,8 +76,46 @@
   logical :: old
 
+  ! help variables for namelist reading
+  integer :: numpoints, parts, readerror
+  integer*2 :: zkind
+  integer :: idate1, itime1, idate2, itime2
+  real :: lon1,lon2,lat1,lat2,z1,z2,mass(nspec)
+  character*40 :: comment
+
+  ! declare namelists
+  namelist /releases_ctrl/ &
+    nspec, &
+    specnum_rel
+
+  namelist /release/ &
+    idate1, itime1, &
+    idate2, itime2, &
+    lon1, lon2, &
+    lat1, lat2, &
+    z1, z2, &
+    zkind, &
+    mass, &
+    parts, &
+    comment
+
+  numpoint=0
+
+  ! presetting namelist releases_ctrl
+  nspec = -1  ! use negative value to determine failed namelist input
+  specnum_rel(1) = 0
+
   !sec, read release to find how many releasepoints should be allocated
-
   open(unitreleases,file=path(1)(1:length(1))//'RELEASES',status='old', &
        err=999)
+
+  ! check if namelist input provided
+  read(unitreleases,releases_ctrl,iostat=readerror)
+
+  if ((readerror.ne.0).or.(nspec.lt.0)) then
+
+  ! no namelist format, close file and allow reopening in old format
+  rewind(unitreleases)
+
+  readerror=1 ! indicates old format
 
   ! Check the format of the RELEASES file (either in free format,
@@ -93,19 +134,15 @@
   rewind(unitreleases)
 
-
   ! Skip first 11 lines (file header)
   !**********************************
 
   call skplin(11,unitreleases)
-
-
   read(unitreleases,*,err=998) nspec
   if (old) call skplin(2,unitreleases)
   do i=1,nspec
-    read(unitreleases,*,err=998) specnum_rel
+    read(unitreleases,*,err=998) specnum_rel(1)
     if (old) call skplin(2,unitreleases)
   end do
 
-  numpoint=0
 100   numpoint=numpoint+1
   read(unitreleases,*,end=25)
@@ -142,50 +179,54 @@
 25   numpoint=numpoint-1
 
+  else 
+    readerror=0
+    do while (readerror.eq.0) 
+      idate1=-1
+      read(unitreleases,release,iostat=readerror)
+      if ((idate1.lt.0).or.(readerror.ne.0)) then
+        readerror=1
+      else
+        numpoint=numpoint+1
+      endif
+    end do
+    readerror=0
+  endif ! if namelist input
+
+  rewind(unitreleases)
+
   !allocate memory for numpoint releaspoint
-    allocate(ireleasestart(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ireleaseend(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xpoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xpoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ypoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(ypoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(zpoint1(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(zpoint2(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(kindz(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xmass(numpoint,maxspec) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(rho_rel(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(npart(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-    allocate(xmasssave(numpoint) &
-         ,stat=stat)
-    if (stat.ne.0) write(*,*)'ERROR: could not allocate RELEASPOINT'
-
-   write (*,*) ' Releasepoints allocated: ', numpoint
-
-   do i=1,numpoint
-     xmasssave(i)=0.
-   end do
+  allocate(ireleasestart(numpoint),stat=stat)
+       
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(ireleaseend(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xpoint1(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xpoint2(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(ypoint1(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(ypoint2(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(zpoint1(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(zpoint2(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(kindz(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xmass(numpoint,maxspec),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(rho_rel(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(npart(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+  allocate(xmasssave(numpoint),stat=stat)
+  if (stat.ne.0) write(*,*)'ERROR: could not allocate releasepoint'
+
+  write (*,*) numpoint,' releasepoints allocated.'
+
+  do i=1,numpoint
+    xmasssave(i)=0.
+  end do
 
   !now save the information
@@ -198,23 +239,25 @@
   end do
 
-  rewind(unitreleases)
-
-
+  if (readerror.ne.0) then
   ! Skip first 11 lines (file header)
   !**********************************
 
-  call skplin(11,unitreleases)
-
+    call skplin(11,unitreleases)
 
   ! Assign species-specific parameters needed for physical processes
   !*************************************************************************
 
-  read(unitreleases,*,err=998) nspec
-  if (nspec.gt.maxspec) goto 994
-  if (old) call skplin(2,unitreleases)
+    read(unitreleases,*,err=998) nspec
+    if (nspec.gt.maxspec) goto 994
+  if (old) call skplin(2,unitreleases)
+  endif
   do i=1,nspec
-    read(unitreleases,*,err=998) specnum_rel
-    if (old) call skplin(2,unitreleases)
-    call readspecies(specnum_rel,i)
+    if (readerror.ne.0) then
+      read(unitreleases,*,err=998) specnum_rel(1)
+      if (old) call skplin(2,unitreleases)
+      call readspecies(specnum_rel(1),i)
+    else
+      call readspecies(specnum_rel(i),i)
+    endif
 
   ! For backward runs, only 1 species is allowed
@@ -301,9 +344,37 @@
   ! Read specifications for each release point
   !*******************************************
-
+  numpoints=numpoint
   numpoint=0
   numpartmax=0
   releaserate=0.
 1000   numpoint=numpoint+1
+
+  if (readerror.eq.0) then ! reading namelist format
+
+  if (numpoint.gt.numpoints) goto 250
+  zkind = 1
+  mass = 0
+  parts = 0
+  comment = ' '
+  read(unitreleases,release,iostat=readerror)
+  id1=idate1
+  it1=itime1
+  id2=idate2
+  it2=itime2
+  xpoint1(numpoint)=lon1
+  xpoint2(numpoint)=lon2
+  ypoint1(numpoint)=lat1
+  ypoint2(numpoint)=lat2
+  zpoint1(numpoint)=z1
+  zpoint2(numpoint)=z2
+  kindz(numpoint)=zkind
+  do i=1,nspec
+    xmass(numpoint,i)=mass(i)
+  end do
+  npart(numpoint)=parts
+  compoint(min(1001,numpoint))=comment
+
+  else
+
   read(unitreleases,*,end=250)
   read(unitreleases,*,err=998,end=250) id1,it1
@@ -347,4 +418,5 @@
   endif
 
+  endif ! if namelist format
 
   ! If a release point contains no particles, stop and issue error message
Index: src_parallel/readspecies.f90
===================================================================
--- src_parallel/readspecies.f90	(revision 3eed2e63a5b06d8d7013ec82d89123c0551cc037)
+++ src_parallel/readspecies.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -30,4 +30,7 @@
   !                                                                            *
   !   11 July 1996                                                             *
+  !
+  !   HSO, 13 August 2013
+  !   added optional namelist input
   !                                                                            *
   !*****************************************************************************
@@ -55,79 +58,116 @@
   logical :: spec_found
 
+  character(len=16) :: pspecies
+  real :: pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer
+  real :: pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
+  integer :: readerror
+
+  ! declare namelist
+  namelist /species_params/ &
+   pspecies, pdecay, pweta, pwetb, preldiff, phenry, pf0, pdensity, pdquer, &
+   pdsigma, pdryvel, pweightmolar, pohreact, pspec_ass, pkao
+
+  pspecies=" "
+  pdecay=-999.9
+  pweta=-9.9E-09
+  pwetb=0.0
+  preldiff=-9.9
+  phenry=0.0
+  pf0=0.0
+  pdensity=-9.9E09
+  pdquer=0.0
+  pdsigma=0.0
+  pdryvel=-9.99
+  pohreact=-9.9E-09
+  pspec_ass=-9
+  pkao=-99.99
+  pweightmolar=-789.0 ! read failure indicator value
+
   ! Open the SPECIES file and read species names and properties
   !************************************************************
   specnum(pos_spec)=id_spec
   write(aspecnumb,'(i3.3)') specnum(pos_spec)
-  open(unitspecies,file= &
-       path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old', &
-       err=998)
-  !write(*,*) 'reading SPECIES',specnum(pos_spec)
+  open(unitspecies,file=path(1)(1:length(1))//'SPECIES/SPECIES_'//aspecnumb,status='old',err=998)
+  write(*,*) 'reading SPECIES',specnum(pos_spec)
 
   ASSSPEC=.FALSE.
 
-  do i=1,6
-    read(unitspecies,*)
-  end do
+  ! try namelist input
+  read(unitspecies,species_params,iostat=readerror)
+   
+  if ((pweightmolar.eq.-789.0).or.(readerror.ne.0)) then ! no namelist found
+
+    rewind(unitspecies)
+
+    do i=1,6
+      read(unitspecies,*)
+    end do
 
     read(unitspecies,'(a10)',end=22) species(pos_spec)
-  !  write(*,*) species(pos_spec)
     read(unitspecies,'(f18.1)',end=22) decay(pos_spec)
-  !  write(*,*) decay(pos_spec)
     read(unitspecies,'(e18.1)',end=22) weta(pos_spec)
-  !  write(*,*) weta(pos_spec)
     read(unitspecies,'(f18.2)',end=22) wetb(pos_spec)
-  !  write(*,*) wetb(pos_spec)
     read(unitspecies,'(f18.1)',end=22) reldiff(pos_spec)
-  !  write(*,*) reldiff(pos_spec)
     read(unitspecies,'(e18.1)',end=22) henry(pos_spec)
-  !  write(*,*) henry(pos_spec)
     read(unitspecies,'(f18.1)',end=22) f0(pos_spec)
-  !  write(*,*) f0(pos_spec)
     read(unitspecies,'(e18.1)',end=22) density(pos_spec)
-  !  write(*,*) density(pos_spec)
     read(unitspecies,'(e18.1)',end=22) dquer(pos_spec)
-  !  write(*,*) dquer(pos_spec)
     read(unitspecies,'(e18.1)',end=22) dsigma(pos_spec)
-  !  write(*,*) dsigma(pos_spec)
     read(unitspecies,'(f18.2)',end=22) dryvel(pos_spec)
-  !  write(*,*) dryvel(pos_spec)
     read(unitspecies,'(f18.2)',end=22) weightmolar(pos_spec)
-  !  write(*,*) weightmolar(pos_spec)
     read(unitspecies,'(e18.1)',end=22) ohreact(pos_spec)
-  !  write(*,*) ohreact(pos_spec)
     read(unitspecies,'(i18)',end=22) spec_ass(pos_spec)
-  !  write(*,*) spec_ass(pos_spec)
     read(unitspecies,'(f18.2)',end=22) kao(pos_spec)
-  !       write(*,*) kao(pos_spec)
-    i=pos_spec
-
-    if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
-       if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
-    endif
-
-    if (spec_ass(pos_spec).gt.0) then
-       spec_found=.FALSE.
-       do j=1,pos_spec-1
-          if (spec_ass(pos_spec).eq.specnum(j)) then
-             spec_ass(pos_spec)=j
-             spec_found=.TRUE.
-             ASSSPEC=.TRUE.
-          endif
-       end do
-       if (spec_found.eqv..FALSE.) then
-          goto 997
-       endif
-    endif
-
-    if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
-    if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
-
-    if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
-      write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
-      write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
-      write(*,*) '#### PARTICLE AND GAS.                       ####'
-      write(*,*) '#### SPECIES NUMBER',aspecnumb
-      stop
-    endif
+
+  else
+
+    species(pos_spec)=pspecies
+    decay(pos_spec)=pdecay
+    weta(pos_spec)=pweta
+    wetb(pos_spec)=pwetb
+    reldiff(pos_spec)=preldiff
+    henry(pos_spec)=phenry
+    f0(pos_spec)=pf0
+    density(pos_spec)=pdensity
+    dquer(pos_spec)=pdquer
+    dsigma(pos_spec)=pdsigma
+    dryvel(pos_spec)=pdryvel
+    weightmolar(pos_spec)=pweightmolar
+    ohreact(pos_spec)=pohreact
+    spec_ass(pos_spec)=pspec_ass
+    kao(pos_spec)=pkao
+
+  endif
+
+  i=pos_spec
+
+  if ((weta(pos_spec).gt.0).and.(henry(pos_spec).le.0)) then
+     if (dquer(pos_spec).le.0) goto 996 ! no particle, no henry set
+  endif
+
+  if (spec_ass(pos_spec).gt.0) then
+     spec_found=.FALSE.
+     do j=1,pos_spec-1
+        if (spec_ass(pos_spec).eq.specnum(j)) then
+           spec_ass(pos_spec)=j
+           spec_found=.TRUE.
+           ASSSPEC=.TRUE.
+        endif
+     end do
+     if (spec_found.eqv..FALSE.) then
+        goto 997
+     endif
+  endif
+
+  if (dsigma(i).eq.1.) dsigma(i)=1.0001   ! avoid floating exception
+  if (dsigma(i).eq.0.) dsigma(i)=1.0001   ! avoid floating exception
+
+  if ((reldiff(i).gt.0.).and.(density(i).gt.0.)) then
+    write(*,*) '#### FLEXPART MODEL ERROR! FILE "SPECIES"    ####'
+    write(*,*) '#### IS CORRUPT. SPECIES CANNOT BE BOTH      ####'
+    write(*,*) '#### PARTICLE AND GAS.                       ####'
+    write(*,*) '#### SPECIES NUMBER',aspecnumb
+    stop
+  endif
 20   continue
 
Index: src_parallel/timemanager.f90
===================================================================
--- src_parallel/timemanager.f90	(revision 7c1fd443e6ad66b1a84c706126014b304c08ff6e)
+++ src_parallel/timemanager.f90	(revision ff050cd28a8d909666c5efc0c10afc7463e6add4)
@@ -97,5 +97,5 @@
        idummy_ran3_initialize, idummy_ran3_advance
 #ifdef NETCDF_OUTPUT
-  use nc_output_mod, only: concoutput_ncdf, concoutput_nest_ncdf
+  use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf
 #endif 
 
@@ -324,7 +324,7 @@
       if ((itime.eq.loutend).and.(outnum.gt.0.)) then
         if ((iout.le.3.).or.(iout.eq.5)) then
-          if (lncdfout) then 
+          if (lnetcdfout.eq.1) then 
 #ifdef NETCDF_OUTPUT
-            call concoutput_ncdf(itime,outnum,gridtotalunc, &
+            call concoutput_netcdf(itime,outnum,gridtotalunc, &
                wetgridtotalunc,drygridtotalunc)
 #endif 
@@ -335,7 +335,7 @@
 
           if (nested_output.eq.1) then
-            if (lncdfout) then 
+            if (lnetcdfout.eq.1) then 
 #ifdef NETCDF_OUTPUT
-              call concoutput_nest_ncdf(itime,outnum)
+              call concoutput_nest_netcdf(itime,outnum)
 #endif 
             else 
