source: flexpart.git/src/calcfluxes.f90

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file was 92fab65, checked in by Ignacio Pisso <ip@…>, 4 years ago

add SPDX-License-Identifier to all .f90 files

  • Property mode set to 100644
File size: 5.9 KB
Line 
1! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
2! SPDX-License-Identifier: GPL-3.0-or-later
3
4subroutine calcfluxes(nage,jpart,xold,yold,zold)
5  !                       i     i    i    i    i
6  !*****************************************************************************
7  !                                                                            *
8  !     Calculation of the gross fluxes across horizontal, eastward and        *
9  !     northward facing surfaces. The routine calculates the mass flux        *
10  !     due to the motion of only one particle. The fluxes of subsequent calls *
11  !     to this subroutine are accumulated until the next output is due.       *
12  !     Upon output, flux fields are re-set to zero in subroutine fluxoutput.f.*
13  !                                                                            *
14  !     Author: A. Stohl                                                       *
15  !                                                                            *
16  !     04 April 2000                                                          *
17  !                                                                            *
18  !*****************************************************************************
19  !                                                                            *
20  ! Variables:                                                                 *
21  !                                                                            *
22  ! nage                  Age class of the particle considered                 *
23  ! jpart                 Index of the particle considered                     *
24  ! xold,yold,zold        "Memorized" old positions of the particle            *
25  !                                                                            *
26  !*****************************************************************************
27
28  use flux_mod
29  use outg_mod
30  use par_mod
31  use com_mod
32
33  implicit none
34
35  integer :: jpart,nage,ixave,jyave,kz,kzave,kp
36  integer :: k,k1,k2,ix,ix1,ix2,ixs,jy,jy1,jy2
37  real :: xold,yold,zold,xmean,ymean
38
39
40  ! Determine average positions
41  !****************************
42
43  if ((ioutputforeachrelease.eq.1).and.(mdomainfill.eq.0)) then
44     kp=npoint(jpart)
45  else
46     kp=1
47  endif
48
49  xmean=(xold+xtra1(jpart))/2.
50  ymean=(yold+ytra1(jpart))/2.
51
52  ixave=int((xmean*dx+xoutshift)/dxout)
53  jyave=int((ymean*dy+youtshift)/dyout)
54  do kz=1,numzgrid                ! determine height of cell
55    if (outheight(kz).gt.ztra1(jpart)) goto 16
56  end do
5716   kzave=kz
58
59
60  ! Determine vertical fluxes
61  !**************************
62
63  if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. &
64       (jyave.le.numygrid-1)) then
65    do kz=1,numzgrid                ! determine height of cell
66      if (outheighthalf(kz).gt.zold) goto 11
67    end do
6811   k1=min(numzgrid,kz)
69    do kz=1,numzgrid                ! determine height of cell
70      if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
71    end do
7221   k2=min(numzgrid,kz)
73
74    do k=1,nspec
75      do kz=k1,k2-1
76        flux(5,ixave,jyave,kz,k,kp,nage)= &
77             flux(5,ixave,jyave,kz,k,kp,nage)+ &
78             xmass1(jpart,k)
79      end do
80      do kz=k2,k1-1
81        flux(6,ixave,jyave,kz,k,kp,nage)= &
82             flux(6,ixave,jyave,kz,k,kp,nage)+ &
83             xmass1(jpart,k)
84      end do
85    end do
86  endif
87
88
89  ! Determine west-east fluxes (fluxw) and east-west fluxes (fluxe)
90  !****************************************************************
91
92  if ((kzave.le.numzgrid).and.(jyave.ge.0).and. &
93       (jyave.le.numygrid-1)) then
94
95  ! 1) Particle does not cross domain boundary
96
97    if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
98      ix1=int((xold*dx+xoutshift)/dxout+0.5)
99      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
100      do k=1,nspec
101        do ix=ix1,ix2-1
102          if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
103            flux(1,ix,jyave,kzave,k,kp,nage)= &
104                 flux(1,ix,jyave,kzave,k,kp,nage) &
105                 +xmass1(jpart,k)
106          endif
107        end do
108        do ix=ix2,ix1-1
109          if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
110            flux(2,ix,jyave,kzave,k,kp,nage)= &
111                 flux(2,ix,jyave,kzave,k,kp,nage) &
112                 +xmass1(jpart,k)
113          endif
114        end do
115      end do
116
117  ! 2) Particle crosses domain boundary: use cyclic boundary condition
118  !    and attribute flux to easternmost grid row only (approximation valid
119  !    for relatively slow motions compared to output grid cell size)
120
121    else
122      ixs=int(((real(nxmin1)-1.e5)*dx+xoutshift)/dxout)
123      if ((ixs.ge.0).and.(ixs.le.numxgrid-1)) then
124        if (xold.gt.xtra1(jpart)) then       ! west-east flux
125          do k=1,nspec
126            flux(1,ixs,jyave,kzave,k,kp,nage)= &
127                 flux(1,ixs,jyave,kzave,k,kp,nage) &
128                 +xmass1(jpart,k)
129          end do
130        else                                 ! east-west flux
131          do k=1,nspec
132            flux(2,ixs,jyave,kzave,k,kp,nage)= &
133                 flux(2,ixs,jyave,kzave,k,kp,nage) &
134                 +xmass1(jpart,k)
135          end do
136        endif
137      endif
138    endif
139  endif
140
141
142  ! Determine south-north fluxes (fluxs) and north-south fluxes (fluxn)
143  !********************************************************************
144
145  if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
146       (ixave.le.numxgrid-1)) then
147    jy1=int((yold*dy+youtshift)/dyout+0.5)
148    jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
149
150    do k=1,nspec
151      do jy=jy1,jy2-1
152        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
153          flux(3,ixave,jy,kzave,k,kp,nage)= &
154               flux(3,ixave,jy,kzave,k,kp,nage) &
155               +xmass1(jpart,k)
156        endif
157      end do
158      do jy=jy2,jy1-1
159        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
160          flux(4,ixave,jy,kzave,k,kp,nage)= &
161               flux(4,ixave,jy,kzave,k,kp,nage) &
162               +xmass1(jpart,k)
163        endif
164      end do
165    end do
166  endif
167
168end subroutine calcfluxes
169
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG