Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/FLEXPART.f90	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
@@ -45,4 +45,9 @@
   use conv_mod
 
+#ifdef NETCDF_OUTPUT
+  use netcdf_output_mod, only: writeheader_netcdf
+#endif
+
+
   implicit none
 
@@ -303,4 +308,22 @@
   !******************************************************************
 
+  if (lnetcdfout.eq.1) then 
+#ifdef NETCDF_OUTPUT
+    call writeheader_netcdf(lnest = .false.)
+#endif
+  else 
+    call writeheader
+  end if
+
+  if (nested_output.eq.1) then
+    if (lnetcdfout.eq.1) then
+#ifdef NETCDF_OUTPUT
+      call writeheader_netcdf(lnest = .true.)
+#endif
+    else
+      call writeheader_nest
+    endif
+  endif
+
   if (verbosity.gt.0) then
     print*,'call writeheader'
@@ -315,5 +338,7 @@
   if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
 
-  !open(unitdates,file=path(2)(1:length(2))//'dates')
+  if (lnetcdfout.ne.1) then
+    open(unitdates,file=path(2)(1:length(2))//'dates')
+  end if
 
   if (verbosity.gt.0) then
@@ -359,9 +384,9 @@
   if (verbosity.gt.0) then
      if (verbosity.gt.1) then   
-       CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
-       write(*,*) 'SYSTEM_CLOCK',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
+       call system_clock(count_clock, count_rate, count_max)
+       write(*,*) 'System clock',(count_clock - count_clock0)/real(count_rate) !, count_rate, count_max
      endif
      if (info_flag.eq.1) then
-       print*, 'info only mode (stop)'    
+       print*, 'Info only mode (stop)'    
        stop
      endif
@@ -373,3 +398,10 @@
   write(*,*) 'CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!'
 
+  ! output wall time
+  if (verbosity .gt. 0) then
+    call system_clock(count_clock,count_rate)
+    tins=(count_clock - count_clock0)/real(count_rate)
+    print*,'Wall time ',tins,'s, ',tins/60,'min, ',tins/3600,'h.'
+  endif
+
 end program flexpart
Index: src/com_mod.f90
===================================================================
--- src/com_mod.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/com_mod.f90	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
@@ -113,4 +113,6 @@
   ! lagespectra  1 if age spectra calculation switched on, 2 if not
 
+  integer :: lnetcdfout
+  ! lnetcdfout   1 for netcdf grid output, 0 if not. Set in COMMAND (namelist input)
 
   integer :: nageclass,lage(maxageclass)
@@ -687,4 +689,5 @@
   integer :: info_flag=0
   integer :: count_clock, count_clock0,  count_rate, count_max
+  real    :: tins
   logical :: nmlout=.true.
    
Index: src/makefile
===================================================================
--- src/makefile	(revision 4fbe7a5609fe43bdc3940768d8f4d949237cc44a)
+++ src/makefile	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
@@ -9,5 +9,4 @@
 #FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
 LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
-
 
 MODOBJS = \
Index: src/makefile.netcdf
===================================================================
--- src/makefile.netcdf	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
+++ src/makefile.netcdf	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
@@ -0,0 +1,95 @@
+SHELL = /bin/bash
+MAIN = FP_ecmwf_gfortran_netcdf
+
+FC       = gfortran
+INCPATH  = /xnilu_wrk/flex_wrk/bin64/grib_api/include
+LIBPATH1 = /xnilu_wrk/flex_wrk/bin64/grib_api/lib
+LIBPATH2 =   /usr/lib/x86_64-linux-gnu/
+NCLIBPATH = /usr/local/lib
+FFLAGS   =   -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+# netcdf output is included by using the complier switch -DNETCDF_OUTPUT and the C preprocessor -cpp
+FFLAGS   =   -O2 -cpp -DNETCDF_OUTPUT -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+#FFLAGS   =   -fbounds-check -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -I$(INCPATH)
+#LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -lgrib_api_f90 -lgrib_api -lm -ljasper
+LDFLAGS  = $(FFLAGS) -L$(LIBPATH2) -L$(LIBPATH1) -L$(NCLIBPATH) -lnetcdf -lnetcdff -lhdf5 -lgrib_api_f90 -lgrib_api -lm -ljasper
+
+MODOBJS = \
+par_mod.o    com_mod.o \
+conv_mod.o              hanna_mod.o \
+interpol_mod.o          cmapf_mod.o \
+unc_mod.o               oh_mod.o \
+xmass_mod.o             flux_mod.o \
+point_mod.o             outg_mod.o \
+netcdf_output_mod.o
+
+OBJECTS = \
+writeheader.o  writeheader_txt.o   writeheader_surf.o       assignland.o\
+calcpar.o               part0.o \
+caldate.o               partdep.o \
+coordtrafo.o            psih.o \
+raerod.o \
+drydepokernel.o         random.o \
+erf.o                   readavailable.o \
+ew.o                    readcommand.o \
+advance.o               readdepo.o \
+releaseparticles.o      psim.o     \
+FLEXPART.o              readlanduse.o \
+getfields.o             init_domainfill.o\
+interpol_wind.o         readoutgrid.o \
+interpol_all.o          readpaths.o \
+getrb.o                 readreceptors.o \
+getrc.o                 readreleases.o \
+getvdep.o               readspecies.o \
+interpol_misslev.o      readwind.o \
+conccalc.o              richardson.o \
+concoutput.o  concoutput_surf.o          scalev.o \
+pbl_profile.o           readOHfield.o\
+juldate.o               timemanager.o \
+interpol_vdep.o         interpol_rain.o \
+verttransform.o         partoutput.o \
+hanna.o                 wetdepokernel.o \
+mean.o                  wetdepo.o \
+hanna_short.o           windalign.o \
+obukhov.o               gridcheck.o \
+hanna1.o                initialize.o \
+                        gridcheck_nests.o \
+readwind_nests.o        calcpar_nests.o \
+verttransform_nests.o   interpol_all_nests.o \
+interpol_wind_nests.o   interpol_misslev_nests.o \
+interpol_vdep_nests.o   interpol_rain_nests.o \
+getvdep_nests.o \
+readageclasses.o        readpartpositions.o \
+calcfluxes.o            fluxoutput.o \
+qvsat.o                 skplin.o \
+convmix.o               calcmatrix.o \
+convect43c.o               redist.o \
+sort2.o                 distance.o \
+centerofmass.o          plumetraj.o \
+openouttraj.o           calcpv.o \
+calcpv_nests.o          distance2.o \
+clustering.o            interpol_wind_short.o \
+interpol_wind_short_nests.o shift_field_0.o \
+shift_field.o           outgrid_init.o \
+openreceptors.o         boundcond_domainfill.o\
+partoutput_short.o      readoutgrid_nest.o \
+outgrid_init_nest.o     writeheader_nest.o writeheader_nest_surf.o \
+concoutput_nest.o concoutput_surf_nest.o      wetdepokernel_nest.o \
+drydepokernel_nest.o    zenithangle.o \
+ohreaction.o            getvdep_nests.o \
+initial_cond_calc.o     initial_cond_output.o \
+dynamic_viscosity.o     get_settling.o
+
+
+$(MAIN): $(MODOBJS) $(OBJECTS)
+	$(FC) *.o -o $(MAIN) $(LDFLAGS)
+
+$(OBJECTS): $(MODOBJS)
+
+%.o: %.f90
+	$(FC) -c $(FFLAGS) $<
+
+clean:
+	rm *.o *.mod
+
+cleanall:
+	rm *.o *.mod $(MAIN)
Index: src/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/readcommand.f90	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
@@ -109,4 +109,5 @@
   nested_output, &
   linit_cond, &
+  lnetcdfout, &
   surf_only    
 
@@ -138,4 +139,5 @@
   nested_output=0
   linit_cond=0
+  lnetcdfout=0
   surf_only=0 
 
@@ -344,4 +346,15 @@
   endif
 
+!  check for netcdf output switch (use for non-namelist input only!)
+  if (iout.ge.8) then
+     lnetcdfout = 1
+     iout = iout - 8
+#ifndef NETCDF_OUTPUT
+     print*,'ERROR: netcdf output not activated during compile time but used in COMMAND file!'
+     print*,'Please recompile with netcdf library or use standard output format.'
+     stop
+#endif
+  endif
+
   ! Check whether a valid option for gridded model output has been chosen
   !**********************************************************************
@@ -349,5 +362,7 @@
   if ((iout.lt.1).or.(iout.gt.5)) then
     write(*,*) ' #### FLEXPART MODEL ERROR! FILE COMMAND:     #### '
-    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5!          #### '
+    write(*,*) ' #### IOUT MUST BE 1, 2, 3, 4, OR 5 FOR       #### '
+    write(*,*) ' #### STANDARD FLEXPART OUTPUT OR  9 - 13     #### '
+    write(*,*) ' #### FOR NETCDF OUTPUT                       #### '
     stop
   endif
@@ -361,5 +376,4 @@
     stop
   endif
-
 
 
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/timemanager.f90	(revision 37017c8f5d1c2f1329bd76c1b7668bf0d4e086cd)
@@ -43,12 +43,14 @@
   !     call convection BEFORE new fields are read in BWD mode                 *
   !  Changes Caroline Forster, Feb 2005                                        *
-  !new interface between flexpart and convection scheme                        *
-  !Emanuel's latest subroutine convect43c.f is used                            *
+  !   new interface between flexpart and convection scheme                     *
+  !   Emanuel's latest subroutine convect43c.f is used                         *
+  !  Changes Stefan Henne, Harald Sodemann, 2013-2014                          *
+  !   added netcdf output code                                                 *
   !*****************************************************************************
   !                                                                            *
   ! Variables:                                                                 *
-  ! DEP                .true. if either wet or dry deposition is switched on   *
+  ! dep                .true. if either wet or dry deposition is switched on   *
   ! decay(maxspec) [1/s] decay constant for radioactive decay                  *
-  ! DRYDEP             .true. if dry deposition is switched on                 *
+  ! drydep             .true. if dry deposition is switched on                 *
   ! ideltas [s]        modelling period                                        *
   ! itime [s]          actual temporal position of calculation                 *
@@ -70,5 +72,5 @@
   ! prob               probability of absorption at ground due to dry          *
   !                    deposition                                              *
-  ! WETDEP             .true. if wet deposition is switched on                 *
+  ! wetdep             .true. if wet deposition is switched on                 *
   ! weight             weight for each concentration sample (1/2 or 1)         *
   ! uap(maxpart),ucp(maxpart),uzp(maxpart) = random velocities due to          *
@@ -92,4 +94,7 @@
   use par_mod
   use com_mod
+#ifdef NETCDF_OUTPUT
+  use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
+#endif 
 
   implicit none
@@ -129,4 +134,5 @@
 
 
+  itime=0
   !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
   write(*,46) float(itime)/3600,itime,numpart
@@ -141,5 +147,4 @@
   do itime=0,ideltas,lsynctime
 
-
   ! Computation of wet deposition, OH reaction and mass transfer
   ! between two species every lsynctime seconds
@@ -152,5 +157,5 @@
   !********************************************************************
 
-    if (WETDEP .and. itime .ne. 0 .and. numpart .gt. 0) then
+    if (wetdep .and. itime .ne. 0 .and. numpart .gt. 0) then
         if (verbosity.gt.0) then
            write (*,*) 'timemanager> call wetdepo'
@@ -159,8 +164,8 @@
     endif
 
-    if (OHREA .and. itime .ne. 0 .and. numpart .gt. 0) &
+    if (ohrea .and. itime .ne. 0 .and. numpart .gt. 0) &
          call ohreaction(itime,lsynctime,loutnext)
 
-    if (ASSSPEC .and. itime .ne. 0 .and. numpart .gt. 0) then
+    if (assspec .and. itime .ne. 0 .and. numpart .gt. 0) then
        stop 'associated species not yet implemented!'
   !     call transferspec(itime,lsynctime,loutnext)
@@ -240,5 +245,5 @@
   !***********************************************************************
 
-    if (DEP.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
+    if (dep.and.(itime.eq.loutnext).and.(ldirect.gt.0)) then
       do ks=1,nspec
       do kp=1,maxpointspec_act
@@ -350,23 +355,48 @@
         if ((iout.le.3.).or.(iout.eq.5)) then
           if (surf_only.ne.1) then 
-          call concoutput(itime,outnum,gridtotalunc, &
-               wetgridtotalunc,drygridtotalunc)
+            if (lnetcdfout.eq.1) then 
+#ifdef NETCDF_OUTPUT
+              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif 
+            else 
+              call concoutput(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+            endif
           else  
-  if (verbosity.eq.1) then
-     print*,'call concoutput_surf '
-     CALL SYSTEM_CLOCK(count_clock)
-     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
-          call concoutput_surf(itime,outnum,gridtotalunc, &
-               wetgridtotalunc,drygridtotalunc)
-  if (verbosity.eq.1) then
-     print*,'called concoutput_surf '
-     CALL SYSTEM_CLOCK(count_clock)
-     WRITE(*,*) 'SYSTEM_CLOCK',count_clock - count_clock0   
-  endif
+            if (verbosity.eq.1) then
+             print*,'call concoutput_surf '
+             call system_clock(count_clock)
+             write(*,*) 'system clock',count_clock - count_clock0   
+            endif
+            if (lnetcdfout.eq.1) then 
+#ifdef NETCDF_OUTPUT
+              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+#endif 
+            else
+              call concoutput_surf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+              if (verbosity.eq.1) then
+                print*,'called concoutput_surf '
+                call system_clock(count_clock)
+                write(*,*) 'system clock',count_clock - count_clock0   
+              endif
+            endif
           endif
 
-          if ((nested_output.eq.1).and.(surf_only.ne.1)) call concoutput_nest(itime,outnum)
-          if ((nested_output.eq.1).and.(surf_only.eq.1)) call concoutput_surf_nest(itime,outnum)
+          if (nested_output .eq. 1) then
+            if (lnetcdfout.eq.0) then
+              if (surf_only.ne.1) then
+                call concoutput_nest(itime,outnum)
+              else 
+                call concoutput_surf_nest(itime,outnum)
+              endif
+            else
+#ifdef NETCDF_OUTPUT
+              if (surf_only.ne.1) then
+                call concoutput_nest_netcdf(itime,outnum)
+              else 
+                call concoutput_surf_nest_netcdf(itime,outnum)
+              endif
+#endif
+            endif
+          endif
           outnum=0.
         endif
@@ -474,6 +504,6 @@
   !****************************
 
-        xold=xtra1(j)
-        yold=ytra1(j)
+        xold=real(xtra1(j))
+        yold=real(ytra1(j))
         zold=ztra1(j)
 
@@ -515,5 +545,5 @@
             endif
 
-            if (DRYDEPSPEC(ks)) then        ! dry deposition
+            if (drydepspec(ks)) then        ! dry deposition
               drydeposit(ks)=xmass1(j,ks)*prob(ks)*decfact
               xmass1(j,ks)=xmass1(j,ks)*(1.-prob(ks))*decfact
@@ -526,5 +556,4 @@
             endif
 
-
             if (mdomainfill.eq.0) then
               if (xmass(npoint(j),ks).gt.0.) &
@@ -538,14 +567,16 @@
           if (xmassfract.lt.0.0001) then   ! terminate all particles carrying less mass
             itra1(j)=-999999999
+            if (verbosity.gt.0) then
+              print*,'terminated particle ',j,' for small mass'
+            endif
           endif
 
   !        Sabine Eckhardt, June 2008
   !        don't create depofield for backward runs
-          if (DRYDEP.AND.(ldirect.eq.1)) then
-            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)), &
-                 real(ytra1(j)),nage,kp)
-            if (nested_output.eq.1) call drydepokernel_nest( &
-                 nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)), &
-                 nage,kp)
+          if (drydep.AND.(ldirect.eq.1)) then
+            call drydepokernel(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
+            if (nested_output.eq.1) then
+              call drydepokernel_nest(nclass(j),drydeposit,real(xtra1(j)),real(ytra1(j)),nage,kp)
+            endif
           endif
 
@@ -554,7 +585,9 @@
 
           if (abs(itra1(j)-itramem(j)).ge.lage(nageclass)) then
-            if (linit_cond.ge.1) &
-                 call initial_cond_calc(itime+lsynctime,j)
+            if (linit_cond.ge.1) call initial_cond_calc(itime+lsynctime,j)
             itra1(j)=-999999999
+            if (verbosity.gt.0) then
+              print*,'terminated particle ',j,' for age'
+            endif
           endif
         endif
@@ -584,11 +617,11 @@
 
   if (iflux.eq.1) then
-      deallocate(flux)
+    deallocate(flux)
   endif
-  if (OHREA.eqv..TRUE.) then
-      deallocate(OH_field,OH_field_height)
+  if (ohrea.eqv..TRUE.) then
+    deallocate(OH_field,OH_field_height)
   endif
   if (ldirect.gt.0) then
-  deallocate(drygridunc,wetgridunc)
+    deallocate(drygridunc,wetgridunc)
   endif
   deallocate(gridunc)
@@ -597,11 +630,11 @@
   deallocate(xmasssave)
   if (nested_output.eq.1) then
-     deallocate(orooutn, arean, volumen)
-     if (ldirect.gt.0) then
-     deallocate(griduncn,drygriduncn,wetgriduncn)
-     endif
+    deallocate(orooutn, arean, volumen)
+    if (ldirect.gt.0) then
+      deallocate(griduncn,drygriduncn,wetgriduncn)
+    endif
   endif
   deallocate(outheight,outheighthalf)
-  deallocate(oroout, area, volume)
+  deallocate(oroout,area,volume)
 
 end subroutine timemanager
