Index: src/fluxoutput.f90
===================================================================
--- src/fluxoutput.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
+++ src/fluxoutput.f90	(revision e200b7adb07d4dfa7bf42f7cb43e0c85b4dbf894)
@@ -0,0 +1,324 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! 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/>.   *
+!**********************************************************************
+
+subroutine fluxoutput(itime)
+  !                        i
+  !*****************************************************************************
+  !                                                                            *
+  !     Output of the gridded fluxes.                                          *
+  !     Eastward, westward, northward, southward, upward and downward gross    *
+  !     fluxes are written to output file in either sparse matrix or grid dump *
+  !     format, whichever is more efficient.                                   *
+  !                                                                            *
+  !     Author: A. Stohl                                                       *
+  !                                                                            *
+  !     04 April 2000                                                          *
+  !                                                                            *
+  !*****************************************************************************
+  !                                                                            *
+  ! Variables:                                                                 *
+  ! ncellse         number of cells with non-zero values for eastward fluxes   *
+  ! sparsee         .true. if in sparse matrix format, else .false.            *
+  !                                                                            *
+  !*****************************************************************************
+
+  use flux_mod
+  use outg_mod
+  use par_mod
+  use com_mod
+
+  implicit none
+
+  real(kind=dp) :: jul
+  integer :: itime,ix,jy,kz,k,nage,jjjjmmdd,ihmmss,kp,i
+  integer :: ncellse(maxspec,maxageclass),ncellsw(maxspec,maxageclass)
+  integer :: ncellss(maxspec,maxageclass),ncellsn(maxspec,maxageclass)
+  integer :: ncellsu(maxspec,maxageclass),ncellsd(maxspec,maxageclass)
+  logical :: sparsee(maxspec,maxageclass),sparsew(maxspec,maxageclass)
+  logical :: sparses(maxspec,maxageclass),sparsen(maxspec,maxageclass)
+  logical :: sparseu(maxspec,maxageclass),sparsed(maxspec,maxageclass)
+  character :: adate*8,atime*6
+
+
+  ! Determine current calendar date, needed for the file name
+  !**********************************************************
+
+  jul=bdate+real(itime,kind=dp)/86400._dp
+  call caldate(jul,jjjjmmdd,ihmmss)
+  write(adate,'(i8.8)') jjjjmmdd
+  write(atime,'(i6.6)') ihmmss
+
+
+  open(unitflux,file=path(2)(1:length(2))//'grid_flux_'//adate// &
+       atime,form='unformatted')
+
+  !**************************************************************
+  ! Check, whether output of full grid or sparse matrix format is
+  ! more efficient in terms of storage space. This is checked for
+  ! every species and for every age class
+  !**************************************************************
+
+  do k=1,nspec
+    do nage=1,nageclass
+      ncellse(k,nage)=0
+      ncellsw(k,nage)=0
+      ncellsn(k,nage)=0
+      ncellss(k,nage)=0
+      ncellsu(k,nage)=0
+      ncellsd(k,nage)=0
+    end do
+  end do
+
+  do k=1,nspec
+  do kp=1,maxpointspec_act
+    do nage=1,nageclass
+      do jy=0,numygrid-1
+        do ix=0,numxgrid-1
+          do kz=1,numzgrid
+            if (flux(2,ix,jy,kz,k,kp,nage).gt.0) ncellse(k,nage)= &
+                 ncellse(k,nage)+1
+            if (flux(1,ix,jy,kz,k,kp,nage).gt.0) ncellsw(k,nage)= &
+                 ncellsw(k,nage)+1
+            if (flux(4,ix,jy,kz,k,kp,nage).gt.0) ncellsn(k,nage)= &
+                 ncellsn(k,nage)+1
+            if (flux(3,ix,jy,kz,k,kp,nage).gt.0) ncellss(k,nage)= &
+                 ncellss(k,nage)+1
+            if (flux(5,ix,jy,kz,k,kp,nage).gt.0) ncellsu(k,nage)= &
+                 ncellsu(k,nage)+1
+            if (flux(6,ix,jy,kz,k,kp,nage).gt.0) ncellsd(k,nage)= &
+                 ncellsd(k,nage)+1
+          end do
+        end do
+      end do
+    end do
+  end do
+  end do
+
+  ! Output in sparse matrix format more efficient, if less than
+  ! 2/5 of all cells contains concentrations>0
+  !************************************************************
+
+  do k=1,nspec
+    do nage=1,nageclass
+      if (4*ncellse(k,nage).lt.numxgrid*numygrid*numzgrid) then
+        sparsee(k,nage)=.true.
+      else
+        sparsee(k,nage)=.false.
+      endif
+      if (4*ncellsw(k,nage).lt.numxgrid*numygrid*numzgrid) then
+        sparsew(k,nage)=.true.
+      else
+        sparsew(k,nage)=.false.
+      endif
+      if (4*ncellsn(k,nage).lt.numxgrid*numygrid*numzgrid) then
+        sparsen(k,nage)=.true.
+      else
+        sparsen(k,nage)=.false.
+      endif
+      if (4*ncellss(k,nage).lt.numxgrid*numygrid*numzgrid) then
+        sparses(k,nage)=.true.
+      else
+        sparses(k,nage)=.false.
+      endif
+      if (4*ncellsu(k,nage).lt.numxgrid*numygrid*numzgrid) then
+        sparseu(k,nage)=.true.
+      else
+        sparseu(k,nage)=.false.
+      endif
+      if (4*ncellsd(k,nage).lt.numxgrid*numygrid*numzgrid) then
+        sparsed(k,nage)=.true.
+      else
+        sparsed(k,nage)=.false.
+      endif
+    end do
+  end do
+
+
+
+  ! Flux output: divide by area and time to get flux in ng/m2/s
+  !************************************************************
+
+  write(unitflux) itime
+  do k=1,nspec
+  do kp=1,maxpointspec_act
+    do nage=1,nageclass
+
+      if (sparsee(k,nage)) then
+        write(unitflux) 1
+        do kz=1,numzgrid
+          do jy=0,numygrid-1
+            do ix=0,numxgrid-1
+              if (flux(2,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
+                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
+                   flux(2,ix,jy,kz,k,kp,nage)/areaeast(ix,jy,kz)/outstep
+            end do
+          end do
+        end do
+        write(unitflux) -999,999.
+      else
+        write(unitflux) 2
+        do kz=1,numzgrid
+          do ix=0,numxgrid-1
+            write(unitflux) (1.e12*flux(2,ix,jy,kz,k,kp,nage)/ &
+                 areaeast(ix,jy,kz)/outstep,jy=0,numygrid-1)
+          end do
+        end do
+      endif
+
+      if (sparsew(k,nage)) then
+        write(unitflux) 1
+        do kz=1,numzgrid
+          do jy=0,numygrid-1
+            do ix=0,numxgrid-1
+              if (flux(1,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
+                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
+                   flux(1,ix,jy,kz,k,kp,nage)/areaeast(ix,jy,kz)/outstep
+            end do
+          end do
+        end do
+        write(unitflux) -999,999.
+      else
+        write(unitflux) 2
+        do kz=1,numzgrid
+          do ix=0,numxgrid-1
+            write(unitflux) (1.e12*flux(1,ix,jy,kz,k,kp,nage)/ &
+                 areaeast(ix,jy,kz)/outstep,jy=0,numygrid-1)
+          end do
+        end do
+      endif
+
+      if (sparses(k,nage)) then
+        write(unitflux) 1
+        do kz=1,numzgrid
+          do jy=0,numygrid-1
+            do ix=0,numxgrid-1
+              if (flux(3,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
+                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
+                   flux(3,ix,jy,kz,k,kp,nage)/areanorth(ix,jy,kz)/outstep
+            end do
+          end do
+        end do
+        write(unitflux) -999,999.
+      else
+        write(unitflux) 2
+        do kz=1,numzgrid
+          do ix=0,numxgrid-1
+            write(unitflux) (1.e12*flux(3,ix,jy,kz,k,kp,nage)/ &
+                 areanorth(ix,jy,kz)/outstep,jy=0,numygrid-1)
+          end do
+        end do
+      endif
+
+      if (sparsen(k,nage)) then
+        write(unitflux) 1
+        do kz=1,numzgrid
+          do jy=0,numygrid-1
+            do ix=0,numxgrid-1 ! north
+              if (flux(4,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
+                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
+                   flux(4,ix,jy,kz,k,kp,nage)/areanorth(ix,jy,kz)/outstep
+            end do
+          end do
+        end do
+        write(unitflux) -999,999.
+      else
+        write(unitflux) 2
+        do kz=1,numzgrid
+          do ix=0,numxgrid-1
+            write(unitflux) (1.e12*flux(4,ix,jy,kz,k,kp,nage)/ &
+                 areanorth(ix,jy,kz)/outstep,jy=0,numygrid-1)
+          end do
+        end do
+      endif
+
+      if (sparseu(k,nage)) then
+        write(unitflux) 1
+        do kz=1,numzgrid
+          do jy=0,numygrid-1
+            do ix=0,numxgrid-1
+              if (flux(5,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
+                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
+                   flux(5,ix,jy,kz,k,kp,nage)/area(ix,jy)/outstep
+            end do
+          end do
+        end do
+        write(unitflux) -999,999.
+      else
+        write(unitflux) 2
+        do kz=1,numzgrid
+          do ix=0,numxgrid-1
+            write(unitflux) (1.e12*flux(5,ix,jy,kz,k,kp,nage)/ &
+                 area(ix,jy)/outstep,jy=0,numygrid-1)
+          end do
+        end do
+      endif
+
+      if (sparsed(k,nage)) then
+        write(unitflux) 1
+        do kz=1,numzgrid
+          do jy=0,numygrid-1
+            do ix=0,numxgrid-1
+              if (flux(6,ix,jy,kz,k,kp,nage).gt.0.) write(unitflux) &
+                   ix+jy*numxgrid+kz*numxgrid*numygrid,1.e12* &
+                   flux(6,ix,jy,kz,k,kp,nage)/area(ix,jy)/outstep
+            end do
+          end do
+        end do
+        write(unitflux) -999,999.
+      else
+        write(unitflux) 2
+        do kz=1,numzgrid
+          do ix=0,numxgrid-1
+            write(unitflux) (1.e12*flux(6,ix,jy,kz,k,kp,nage)/ &
+                 area(ix,jy)/outstep,jy=0,numygrid-1)
+          end do
+        end do
+      endif
+
+    end do
+  end do
+  end do
+
+
+  close(unitflux)
+
+
+  ! Reinitialization of grid
+  !*************************
+
+  do k=1,nspec
+  do kp=1,maxpointspec_act
+    do jy=0,numygrid-1
+      do ix=0,numxgrid-1
+          do kz=1,numzgrid
+            do nage=1,nageclass
+              do i=1,6
+                flux(i,ix,jy,kz,k,kp,nage)=0.
+              end do
+            end do
+          end do
+      end do
+    end do
+  end do
+  end do
+
+
+end subroutine fluxoutput
