Index: src/FLEXPART.f90
===================================================================
--- src/FLEXPART.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/FLEXPART.f90	(revision 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -45,4 +45,7 @@
   use conv_mod
 
+  use netcdf_output_mod, only: writeheader_netcdf
+
+
   implicit none
 
@@ -303,4 +306,18 @@
   !******************************************************************
 
+  if (lnetcdfout.eq.1) then 
+    call writeheader_netcdf(lnest = .false.)
+  else 
+    call writeheader
+  end if
+
+  if (nested_output.eq.1) then
+    if (lnetcdfout.eq.1) then
+      call writeheader_netcdf(lnest = .true.)
+    else
+      call writeheader_nest
+    endif
+  endif
+
   if (verbosity.gt.0) then
     print*,'call writeheader'
@@ -315,6 +332,4 @@
   if (nested_output.ne.1.and.surf_only.eq.1) call writeheader_surf
 
-  !open(unitdates,file=path(2)(1:length(2))//'dates')
-
   if (verbosity.gt.0) then
     print*,'call openreceptors'
@@ -359,9 +374,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 +388,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 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -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/readcommand.f90
===================================================================
--- src/readcommand.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/readcommand.f90	(revision 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -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,10 @@
   endif
 
+!  check for netcdf output switch (use for non-namelist input only!)
+  if (iout.ge.8) then
+     lnetcdfout = 1
+     iout = iout - 8
+  endif
+
   ! Check whether a valid option for gridded model output has been chosen
   !**********************************************************************
@@ -349,5 +357,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 +371,4 @@
     stop
   endif
-
 
 
Index: src/readreceptors.f90
===================================================================
--- src/readreceptors.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/readreceptors.f90	(revision 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -81,5 +81,5 @@
   ! prepare namelist output if requested
   if (nmlout.eqv..true.) then
-    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',access='append',status='new',err=1000)
+    open(unitreceptorout,file=path(2)(1:length(2))//'RECEPTORS.namelist',access='append',status='unknown',err=1000)
   endif
 
@@ -179,6 +179,6 @@
   stop
 
-1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS"    #### '
-  write(*,*) ' #### CANNOT BE OPENED IN THE DIRECTORY       #### '
+1000 write(*,*) ' #### FLEXPART MODEL ERROR! FILE "RECEPTORS.namelist"    #### '
+  write(*,*) ' #### CANNOT BE OPENED/CREATED IN THE DIRECTORY       #### '
   write(*,'(a)') path(2)(1:length(2))
   stop
Index: src/readreleases.f90
===================================================================
--- src/readreleases.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/readreleases.f90	(revision 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -126,5 +126,5 @@
   ! prepare namelist output if requested
   if (nmlout.eqv..true.) then
-    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='new',err=1000)
+    open(unitreleasesout,file=path(2)(1:length(2))//'RELEASES.namelist',access='append',status='unknown',err=1000)
   endif
 
Index: src/readspecies.f90
===================================================================
--- src/readspecies.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/readspecies.f90	(revision 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -239,18 +239,20 @@
 20   continue
 
+
+  ! Read in daily and day-of-week variation of emissions, if available
+  !*******************************************************************
+  ! HSO: This is not yet implemented as namelist parameters
+  ! default values set to 1
+
+  do j=1,24           ! initialize everything to no variation
+    area_hour(i,j)=1.
+    point_hour(i,j)=1.
+  end do
+  do j=1,7
+    area_dow(i,j)=1.
+    point_dow(i,j)=1.
+  end do
+
   if (readerror.ne.0) then ! text format input
-
-    ! Read in daily and day-of-week variation of emissions, if available
-    !*******************************************************************
-    ! HSO: This is not yet implemented as namelist parameters
-
-    do j=1,24           ! initialize everything to no variation
-      area_hour(i,j)=1.
-      point_hour(i,j)=1.
-    end do
-    do j=1,7
-      area_dow(i,j)=1.
-      point_dow(i,j)=1.
-    end do
 
     read(unitspecies,*,end=22)
@@ -269,5 +271,5 @@
   ! namelist output if requested
   if (nmlout.eqv..true.) then
-    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='new',err=1000)
+    open(unitspecies,file=path(2)(1:length(2))//'SPECIES_'//aspecnumb//'.namelist',access='append',status='unknown',err=1000)
     write(unitspecies,nml=species_params)
     close(unitspecies)
Index: src/timemanager.f90
===================================================================
--- src/timemanager.f90	(revision b4d29ce9942e7817568d557b23ac2d4ccb3bfb37)
+++ src/timemanager.f90	(revision 660bcb7e224bffb410c0aac2ae9d9831df216971)
@@ -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,5 @@
   use par_mod
   use com_mod
+  use netcdf_output_mod, only: concoutput_netcdf, concoutput_nest_netcdf,concoutput_surf_netcdf, concoutput_surf_nest_netcdf
 
   implicit none
@@ -129,4 +132,5 @@
 
 
+  itime=0
   !write(*,45) itime,numpart,gridtotalunc,wetgridtotalunc,drygridtotalunc
   write(*,46) float(itime)/3600,itime,numpart
@@ -141,5 +145,4 @@
   do itime=0,ideltas,lsynctime
 
-
   ! Computation of wet deposition, OH reaction and mass transfer
   ! between two species every lsynctime seconds
@@ -152,5 +155,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 +162,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 +243,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 +353,42 @@
         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 
+              call concoutput_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+            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 
+              call concoutput_surf_netcdf(itime,outnum,gridtotalunc,wetgridtotalunc,drygridtotalunc)
+            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
+              if (surf_only.ne.1) then
+                call concoutput_nest_netcdf(itime,outnum)
+              else 
+                call concoutput_surf_nest_netcdf(itime,outnum)
+              endif
+            endif
+          endif
           outnum=0.
         endif
@@ -474,6 +496,6 @@
   !****************************
 
-        xold=xtra1(j)
-        yold=ytra1(j)
+        xold=real(xtra1(j))
+        yold=real(ytra1(j))
         zold=ztra1(j)
 
@@ -515,5 +537,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 +548,4 @@
             endif
 
-
             if (mdomainfill.eq.0) then
               if (xmass(npoint(j),ks).gt.0.) &
@@ -538,14 +559,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 +577,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 +609,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 +622,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
