! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later module flux_mod ! flux eastward, westward, northward, southward, upward and downward ! fluxes of all species and all ageclasses ! areaeast,areanorth [m2] side areas of each grid cell use outgrid_mod use par_mod use com_mod use windfields_mod implicit none !Moved to outgrid_mod, because of dependencies ! real,allocatable, dimension (:,:,:,:,:,:,:) :: flux !1 fluxw west - east !2 fluxe east - west !3 fluxs south - north !4 fluxn north - south !5 fluxu upward !6 fluxd downward !real,allocatable, dimension (:,:,:) :: areanorth !real,allocatable, dimension (:,:,:) :: areaeast contains subroutine calcfluxes(itime,nage,jpart,xold,yold,zold,thread) ! i i i i i !***************************************************************************** ! * ! Calculation of the gross fluxes across horizontal, eastward and * ! northward facing surfaces. The routine calculates the mass flux * ! due to the motion of only one particle. The fluxes of subsequent calls * ! to this subroutine are accumulated until the next output is due. * ! Upon output, flux fields are re-set to zero in subroutine fluxoutput.f.* ! * ! Author: A. Stohl * ! * ! 04 April 2000 * ! * ! Changes * ! 2021 L. Bakels: OpenMP parallelisation * ! * !***************************************************************************** ! * ! Variables: * ! * ! nage Age class of the particle considered * ! jpart Index of the particle considered * ! xold,yold,zold "Memorized" old positions of the particle * ! * !***************************************************************************** use particle_mod #ifdef ETA use coord_ecmwf_mod #endif implicit none integer, intent(in) :: thread ! for OMP, number of thread integer :: itime,jpart,nage,ixave,jyave,kz,kzave,kp integer :: k,k1,k2,ix,ix1,ix2,ixs,jy,jy1,jy2 real :: xold,yold,zold,xmean,ymean ! Determine average positions !**************************** if ((ioutputforeachrelease.eq.1).and.(mdomainfill.eq.0)) then kp=part(jpart)%npoint else kp=1 endif #ifdef ETA call update_zeta_to_z(itime,jpart) #endif xmean=(xold+real(part(jpart)%xlon))*0.5 ymean=(yold+real(part(jpart)%ylat))*0.5 ixave=int((xmean*dx+xoutshift)/dxout) jyave=int((ymean*dy+youtshift)/dyout) do kz=1,numzgrid ! determine height of cell if (outheight(kz).gt.part(jpart)%z) exit end do kzave=kz ! Determine vertical fluxes !************************** if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. & (jyave.le.numygrid-1)) then do kz=1,numzgrid ! determine height of cell if (outheighthalf(kz).gt.zold) exit end do k1=min(numzgrid,kz) do kz=1,numzgrid ! determine height of cell if (outheighthalf(kz).gt.part(jpart)%z) exit end do k2=min(numzgrid,kz) do k=1,nspec do kz=k1,k2-1 #ifdef _OPENMP flux_omp(5,ixave,jyave,kz,k,kp,nage,thread)= & flux_omp(5,ixave,jyave,kz,k,kp,nage,thread)+ & mass(jpart,k) #else flux(5,ixave,jyave,kz,k,kp,nage)= & flux(5,ixave,jyave,kz,k,kp,nage)+ & mass(jpart,k) #endif end do do kz=k2,k1-1 #ifdef _OPENMP flux_omp(6,ixave,jyave,kz,k,kp,nage,thread)= & flux_omp(6,ixave,jyave,kz,k,kp,nage,thread)+ & mass(jpart,k) #else flux(6,ixave,jyave,kz,k,kp,nage)= & flux(6,ixave,jyave,kz,k,kp,nage)+ & mass(jpart,k) #endif end do end do endif ! Determine west-east fluxes (fluxw) and east-west fluxes (fluxe) !**************************************************************** if ((kzave.le.numzgrid).and.(jyave.ge.0).and. & (jyave.le.numygrid-1)) then ! 1) Particle does not cross domain boundary if (abs(xold-part(jpart)%xlon).lt.real(nx)*0.5) then ix1=int((xold*dx+xoutshift)/dxout+0.5) ix2=int((part(jpart)%xlon*dx+xoutshift)/dxout+0.5) do k=1,nspec do ix=ix1,ix2-1 if ((ix.ge.0).and.(ix.le.numxgrid-1)) then #ifdef _OPENMP flux_omp(1,ix,jyave,kzave,k,kp,nage,thread)= & flux_omp(1,ix,jyave,kzave,k,kp,nage,thread) & +mass(jpart,k) #else flux(1,ix,jyave,kzave,k,kp,nage)= & flux(1,ix,jyave,kzave,k,kp,nage) & +mass(jpart,k) #endif endif end do do ix=ix2,ix1-1 if ((ix.ge.0).and.(ix.le.numxgrid-1)) then #ifdef _OPENMP flux_omp(2,ix,jyave,kzave,k,kp,nage,thread)= & flux_omp(2,ix,jyave,kzave,k,kp,nage,thread) & +mass(jpart,k) #else flux(2,ix,jyave,kzave,k,kp,nage)= & flux(2,ix,jyave,kzave,k,kp,nage) & +mass(jpart,k) #endif endif end do end do ! 2) Particle crosses domain boundary: use cyclic boundary condition ! and attribute flux to easternmost grid row only (approximation valid ! for relatively slow motions compared to output grid cell size) else ixs=int(((real(nxmin1)-1.e5)*dx+xoutshift)/dxout) if ((ixs.ge.0).and.(ixs.le.numxgrid-1)) then if (xold.gt.part(jpart)%xlon) then ! west-east flux do k=1,nspec #ifdef _OPENMP flux_omp(1,ixs,jyave,kzave,k,kp,nage,thread)= & flux_omp(1,ixs,jyave,kzave,k,kp,nage,thread) & +mass(jpart,k) #else flux(1,ixs,jyave,kzave,k,kp,nage)= & flux(1,ixs,jyave,kzave,k,kp,nage) & +mass(jpart,k) #endif end do else ! east-west flux do k=1,nspec #ifdef _OPENMP flux_omp(2,ixs,jyave,kzave,k,kp,nage,thread)= & flux_omp(2,ixs,jyave,kzave,k,kp,nage,thread) & +mass(jpart,k) #else flux(2,ixs,jyave,kzave,k,kp,nage)= & flux(2,ixs,jyave,kzave,k,kp,nage) & +mass(jpart,k) #endif end do endif endif endif endif ! Determine south-north fluxes (fluxs) and north-south fluxes (fluxn) !******************************************************************** if ((kzave.le.numzgrid).and.(ixave.ge.0).and. & (ixave.le.numxgrid-1)) then jy1=int((yold*dy+youtshift)/dyout+0.5) jy2=int((part(jpart)%ylat*dy+youtshift)/dyout+0.5) do k=1,nspec do jy=jy1,jy2-1 if ((jy.ge.0).and.(jy.le.numygrid-1)) then #ifdef _OPENMP flux_omp(3,ixave,jy,kzave,k,kp,nage,thread)= & flux_omp(3,ixave,jy,kzave,k,kp,nage,thread) & +mass(jpart,k) #else flux(3,ixave,jy,kzave,k,kp,nage)= & flux(3,ixave,jy,kzave,k,kp,nage) & +mass(jpart,k) #endif endif end do do jy=jy2,jy1-1 if ((jy.ge.0).and.(jy.le.numygrid-1)) then #ifdef _OPENMP flux_omp(4,ixave,jy,kzave,k,kp,nage,thread)= & flux_omp(4,ixave,jy,kzave,k,kp,nage,thread) & +mass(jpart,k) #else flux(4,ixave,jy,kzave,k,kp,nage)= & flux(4,ixave,jy,kzave,k,kp,nage) & +mass(jpart,k) #endif endif end do end do endif end subroutine calcfluxes 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 date_mod implicit none real(kind=dp) :: jul integer :: itime,ix,jy,kz,k,nage,jjjjmmdd,ihmmss,kp,i integer :: ncellse(maxspec,nageclass),ncellsw(maxspec,nageclass) integer :: ncellss(maxspec,nageclass),ncellsn(maxspec,nageclass) integer :: ncellsu(maxspec,nageclass),ncellsd(maxspec,nageclass) logical :: sparsee(maxspec,nageclass),sparsew(maxspec,nageclass) logical :: sparses(maxspec,nageclass),sparsen(maxspec,nageclass) logical :: sparseu(maxspec,nageclass),sparsed(maxspec,nageclass) 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) write(*,*) 'Flux:', itime, flux(:,1,1,1,1,1,1) ! 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. #ifdef _OPENMP flux_omp(i,ix,jy,kz,k,kp,nage,:)=0. #endif end do end do end do end do end do end do end do end subroutine fluxoutput end module flux_mod