fluxoutput.f90

Go to the documentation of this file.

!**********************************************************************
! 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