! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt ! SPDX-License-Identifier: GPL-3.0-or-later subroutine calcfluxes(nage,jpart,xold,yold,zold) ! 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 * ! * !***************************************************************************** ! * ! Variables: * ! * ! nage Age class of the particle considered * ! jpart Index of the particle considered * ! xold,yold,zold "Memorized" old positions of the particle * ! * !***************************************************************************** use flux_mod use outg_mod use par_mod use com_mod implicit none integer :: 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=npoint(jpart) else kp=1 endif xmean=(xold+xtra1(jpart))/2. ymean=(yold+ytra1(jpart))/2. ixave=int((xmean*dx+xoutshift)/dxout) jyave=int((ymean*dy+youtshift)/dyout) do kz=1,numzgrid ! determine height of cell if (outheight(kz).gt.ztra1(jpart)) goto 16 end do 16 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) goto 11 if (outheight(kz).gt.zold) goto 11 !sec, use upper layer instead of mid layer end do 11 k1=min(numzgrid,kz) do kz=1,numzgrid ! determine height of cell ! if (outheighthalf(kz).gt.ztra1(jpart)) goto 21 if (outheight(kz).gt.ztra1(jpart)) goto 21 end do 21 k2=min(numzgrid,kz) do k=1,nspec do kz=k1,k2-1 flux(5,ixave,jyave,kz,k,kp,nage)= & flux(5,ixave,jyave,kz,k,kp,nage)+ & xmass1(jpart,k) end do do kz=k2,k1-1 flux(6,ixave,jyave,kz,k,kp,nage)= & flux(6,ixave,jyave,kz,k,kp,nage)+ & xmass1(jpart,k) 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-xtra1(jpart)).lt.real(nx)/2.) then ix1=int((xold*dx+xoutshift)/dxout) ! flux throught the western boundary (sec) ix2=int((xtra1(jpart)*dx+xoutshift)/dxout) ! flux throught the western boundary (sec) ! ix1=int((xold*dx+xoutshift)/dxout+0.5) ! ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5) do k=1,nspec do ix=ix1,ix2-1 if ((ix.ge.0).and.(ix.le.numxgrid-1)) then flux(1,ix,jyave,kzave,k,kp,nage)= & flux(1,ix,jyave,kzave,k,kp,nage) & +xmass1(jpart,k) endif end do do ix=ix2,ix1-1 if ((ix.ge.0).and.(ix.le.numxgrid-1)) then flux(2,ix,jyave,kzave,k,kp,nage)= & flux(2,ix,jyave,kzave,k,kp,nage) & +xmass1(jpart,k) 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.xtra1(jpart)) then ! west-east flux do k=1,nspec flux(1,ixs,jyave,kzave,k,kp,nage)= & flux(1,ixs,jyave,kzave,k,kp,nage) & +xmass1(jpart,k) end do else ! east-west flux do k=1,nspec flux(2,ixs,jyave,kzave,k,kp,nage)= & flux(2,ixs,jyave,kzave,k,kp,nage) & +xmass1(jpart,k) 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) ! flux throught the southern boundary (sec) jy2=int((ytra1(jpart)*dy+youtshift)/dyout) ! flux throught the southern boundary (sec) ! jy1=int((yold*dy+youtshift)/dyout+0.5) ! jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5) do k=1,nspec do jy=jy1,jy2-1 if ((jy.ge.0).and.(jy.le.numygrid-1)) then flux(3,ixave,jy,kzave,k,kp,nage)= & flux(3,ixave,jy,kzave,k,kp,nage) & +xmass1(jpart,k) endif end do do jy=jy2,jy1-1 if ((jy.ge.0).and.(jy.le.numygrid-1)) then flux(4,ixave,jy,kzave,k,kp,nage)= & flux(4,ixave,jy,kzave,k,kp,nage) & +xmass1(jpart,k) endif end do end do endif end subroutine calcfluxes