source: branches/flexpart91_hasod/src_parallel/calcfluxes.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 9 years ago

Added parallel version of Flexpart91

File size: 7.6 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine calcfluxes(nage,jpart,xold,yold,zold)
23  !                       i     i    i    i    i
24  !*****************************************************************************
25  !                                                                            *
26  !     Calculation of the gross fluxes across horizontal, eastward and        *
27  !     northward facing surfaces. The routine calculates the mass flux        *
28  !     due to the motion of only one particle. The fluxes of subsequent calls *
29  !     to this subroutine are accumulated until the next output is due.       *
30  !     Upon output, flux fields are re-set to zero in subroutine fluxoutput.f.*
31  !                                                                            *
32  !     Author: A. Stohl                                                       *
33  !                                                                            *
34  !     04 April 2000                                                          *
35  !                                                                            *
36  !*****************************************************************************
37  !                                                                            *
38  ! Variables:                                                                 *
39  !                                                                            *
40  ! nage                  Age class of the particle considered                 *
41  ! jpart                 Index of the particle considered                     *
42  ! xold,yold,zold        "Memorized" old positions of the particle            *
43  !                                                                            *
44  !*****************************************************************************
45
46  use flux_mod, only: flux
47  use outg_mod, only: outheighthalf, outheight
48  use com_mod, only: nx, dx, numzgrid, xmass1, xtra1, ytra1, ztra1, npoint, xoutshift, &
49    youtshift, nxmin1, nymin1, ny, dy, numxgrid, numygrid, nspec, ioutputforeachrelease, &
50    dxout, dyout, mdomainfill
51
52  implicit none
53
54  integer :: jpart,nage,ixave,jyave,kz,kzave,kp
55  integer :: k,k1,k2,ix,ix1,ix2,ixs,jy,jy1,jy2
56  real :: xold,yold,zold,xmean,ymean
57
58
59  ! Determine average positions
60  !****************************
61
62  if ((ioutputforeachrelease.eq.1).and.(mdomainfill.eq.0)) then
63     kp=npoint(jpart)
64  else
65     kp=1
66  endif
67
68  xmean=(xold+real(xtra1(jpart)))/2.
69  ymean=(yold+real(ytra1(jpart)))/2.
70
71  ixave=int((xmean*dx+xoutshift)/dxout)
72  jyave=int((ymean*dy+youtshift)/dyout)
73  do kz=1,numzgrid                ! determine height of cell
74    if (outheight(kz).gt.ztra1(jpart)) goto 16
75  end do
7616   kzave=kz
77
78
79  ! Determine vertical fluxes
80  !**************************
81
82  if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. &
83       (jyave.le.numygrid-1)) then
84    do kz=1,numzgrid                ! determine height of cell
85      if (outheighthalf(kz).gt.zold) goto 11
86    end do
8711   k1=min(numzgrid,kz)
88    do kz=1,numzgrid                ! determine height of cell
89      if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
90    end do
9121   k2=min(numzgrid,kz)
92
93!$OMP CRITICAL
94    do k=1,nspec
95      do kz=k1,k2-1
96        flux(5,ixave,jyave,kz,k,kp,nage)= &
97             flux(5,ixave,jyave,kz,k,kp,nage)+ &
98             xmass1(jpart,k)
99      end do
100      do kz=k2,k1-1
101        flux(6,ixave,jyave,kz,k,kp,nage)= &
102             flux(6,ixave,jyave,kz,k,kp,nage)+ &
103             xmass1(jpart,k)
104      end do
105    end do
106!$OMP END CRITICAL
107  endif
108
109
110  ! Determine west-east fluxes (fluxw) and east-west fluxes (fluxe)
111  !****************************************************************
112
113  if ((kzave.le.numzgrid).and.(jyave.ge.0).and. &
114       (jyave.le.numygrid-1)) then
115
116  ! 1) Particle does not cross domain boundary
117
118    if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
119      ix1=int((xold*dx+xoutshift)/dxout+0.5)
120      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
121!$OMP CRITICAL
122      do k=1,nspec
123        do ix=ix1,ix2-1
124          if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
125            flux(1,ix,jyave,kzave,k,kp,nage)= &
126                 flux(1,ix,jyave,kzave,k,kp,nage) &
127                 +xmass1(jpart,k)
128          endif
129        end do
130        do ix=ix2,ix1-1
131          if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
132            flux(2,ix,jyave,kzave,k,kp,nage)= &
133                 flux(2,ix,jyave,kzave,k,kp,nage) &
134                 +xmass1(jpart,k)
135          endif
136        end do
137      end do
138!$OMP END CRITICAL
139
140  ! 2) Particle crosses domain boundary: use cyclic boundary condition
141  !    and attribute flux to easternmost grid row only (approximation valid
142  !    for relatively slow motions compared to output grid cell size)
143
144    else
145      ixs=int(((real(nxmin1)-1.e5)*dx+xoutshift)/dxout)
146      if ((ixs.ge.0).and.(ixs.le.numxgrid-1)) then
147!$OMP CRITICAL
148        if (xold.gt.xtra1(jpart)) then       ! west-east flux
149          do k=1,nspec
150            flux(1,ixs,jyave,kzave,k,kp,nage)= &
151                 flux(1,ixs,jyave,kzave,k,kp,nage) &
152                 +xmass1(jpart,k)
153          end do
154        else                                 ! east-west flux
155          do k=1,nspec
156            flux(2,ixs,jyave,kzave,k,kp,nage)= &
157                 flux(2,ixs,jyave,kzave,k,kp,nage) &
158                 +xmass1(jpart,k)
159          end do
160        endif
161!$OMP END CRITICAL
162      endif
163    endif
164  endif
165
166
167  ! Determine south-north fluxes (fluxs) and north-south fluxes (fluxn)
168  !********************************************************************
169
170  if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
171       (ixave.le.numxgrid-1)) then
172    jy1=int((yold*dy+youtshift)/dyout+0.5)
173    jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
174
175!$OMP CRITICAL
176    do k=1,nspec
177      do jy=jy1,jy2-1
178        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
179          flux(3,ixave,jy,kzave,k,kp,nage)= &
180               flux(3,ixave,jy,kzave,k,kp,nage) &
181               +xmass1(jpart,k)
182        endif
183      end do
184      do jy=jy2,jy1-1
185        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
186          flux(4,ixave,jy,kzave,k,kp,nage)= &
187               flux(4,ixave,jy,kzave,k,kp,nage) &
188               +xmass1(jpart,k)
189        endif
190      end do
191    end do
192!$OMP END CRITICAL
193  endif
194
195end subroutine calcfluxes
196
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG