source: trunk/src/calcfluxes.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 9 years ago
File size: 7.2 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
47  use outg_mod
48  use par_mod
49  use com_mod
50
51  implicit none
52
53  integer :: jpart,nage,ixave,jyave,kz,kzave,kp
54  integer :: k,k1,k2,ix,ix1,ix2,ixs,jy,jy1,jy2
55  real :: xold,yold,zold,xmean,ymean
56
57
58  ! Determine average positions
59  !****************************
60
61  if ((ioutputforeachrelease.eq.1).and.(mdomainfill.eq.0)) then
62     kp=npoint(jpart)
63  else
64     kp=1
65  endif
66
67  xmean=(xold+xtra1(jpart))/2.
68  ymean=(yold+ytra1(jpart))/2.
69
70  ixave=int((xmean*dx+xoutshift)/dxout)
71  jyave=int((ymean*dy+youtshift)/dyout)
72  do kz=1,numzgrid                ! determine height of cell
73    if (outheight(kz).gt.ztra1(jpart)) goto 16
74  end do
7516   kzave=kz
76
77
78  ! Determine vertical fluxes
79  !**************************
80
81  if ((ixave.ge.0).and.(jyave.ge.0).and.(ixave.le.numxgrid-1).and. &
82       (jyave.le.numygrid-1)) then
83    do kz=1,numzgrid                ! determine height of cell
84      if (outheighthalf(kz).gt.zold) goto 11
85    end do
8611   k1=min(numzgrid,kz)
87    do kz=1,numzgrid                ! determine height of cell
88      if (outheighthalf(kz).gt.ztra1(jpart)) goto 21
89    end do
9021   k2=min(numzgrid,kz)
91
92    do k=1,nspec
93      do kz=k1,k2-1
94        flux(5,ixave,jyave,kz,k,kp,nage)= &
95             flux(5,ixave,jyave,kz,k,kp,nage)+ &
96             xmass1(jpart,k)
97      end do
98      do kz=k2,k1-1
99        flux(6,ixave,jyave,kz,k,kp,nage)= &
100             flux(6,ixave,jyave,kz,k,kp,nage)+ &
101             xmass1(jpart,k)
102      end do
103    end do
104  endif
105
106
107  ! Determine west-east fluxes (fluxw) and east-west fluxes (fluxe)
108  !****************************************************************
109
110  if ((kzave.le.numzgrid).and.(jyave.ge.0).and. &
111       (jyave.le.numygrid-1)) then
112
113  ! 1) Particle does not cross domain boundary
114
115    if (abs(xold-xtra1(jpart)).lt.real(nx)/2.) then
116      ix1=int((xold*dx+xoutshift)/dxout+0.5)
117      ix2=int((xtra1(jpart)*dx+xoutshift)/dxout+0.5)
118      do k=1,nspec
119        do ix=ix1,ix2-1
120          if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
121            flux(1,ix,jyave,kzave,k,kp,nage)= &
122                 flux(1,ix,jyave,kzave,k,kp,nage) &
123                 +xmass1(jpart,k)
124          endif
125        end do
126        do ix=ix2,ix1-1
127          if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
128            flux(2,ix,jyave,kzave,k,kp,nage)= &
129                 flux(2,ix,jyave,kzave,k,kp,nage) &
130                 +xmass1(jpart,k)
131          endif
132        end do
133      end do
134
135  ! 2) Particle crosses domain boundary: use cyclic boundary condition
136  !    and attribute flux to easternmost grid row only (approximation valid
137  !    for relatively slow motions compared to output grid cell size)
138
139    else
140      ixs=int(((real(nxmin1)-1.e5)*dx+xoutshift)/dxout)
141      if ((ixs.ge.0).and.(ixs.le.numxgrid-1)) then
142        if (xold.gt.xtra1(jpart)) then       ! west-east flux
143          do k=1,nspec
144            flux(1,ixs,jyave,kzave,k,kp,nage)= &
145                 flux(1,ixs,jyave,kzave,k,kp,nage) &
146                 +xmass1(jpart,k)
147          end do
148        else                                 ! east-west flux
149          do k=1,nspec
150            flux(2,ixs,jyave,kzave,k,kp,nage)= &
151                 flux(2,ixs,jyave,kzave,k,kp,nage) &
152                 +xmass1(jpart,k)
153          end do
154        endif
155      endif
156    endif
157  endif
158
159
160  ! Determine south-north fluxes (fluxs) and north-south fluxes (fluxn)
161  !********************************************************************
162
163  if ((kzave.le.numzgrid).and.(ixave.ge.0).and. &
164       (ixave.le.numxgrid-1)) then
165    jy1=int((yold*dy+youtshift)/dyout+0.5)
166    jy2=int((ytra1(jpart)*dy+youtshift)/dyout+0.5)
167
168    do k=1,nspec
169      do jy=jy1,jy2-1
170        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
171          flux(3,ixave,jy,kzave,k,kp,nage)= &
172               flux(3,ixave,jy,kzave,k,kp,nage) &
173               +xmass1(jpart,k)
174        endif
175      end do
176      do jy=jy2,jy1-1
177        if ((jy.ge.0).and.(jy.le.numygrid-1)) then
178          flux(4,ixave,jy,kzave,k,kp,nage)= &
179               flux(4,ixave,jy,kzave,k,kp,nage) &
180               +xmass1(jpart,k)
181        endif
182      end do
183    end do
184  endif
185
186end subroutine calcfluxes
187
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG