source: trunk/src/conccalc.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 14.8 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 conccalc(itime,weight)
23  !                      i     i
24  !*****************************************************************************
25  !                                                                            *
26  !     Calculation of the concentrations on a regular grid using volume       *
27  !     sampling                                                               *
28  !                                                                            *
29  !     Author: A. Stohl                                                       *
30  !                                                                            *
31  !     24 May 1996                                                            *
32  !                                                                            *
33  !     April 2000: Update to calculate age spectra                            *
34  !                 Bug fix to avoid negative conc. at the domain boundaries,  *
35  !                 as suggested by Petra Seibert                              *
36  !                                                                            *
37  !     2 July 2002: re-order if-statements in order to optimize CPU time      *
38  !                                                                            *
39  !                                                                            *
40  !*****************************************************************************
41  !                                                                            *
42  ! Variables:                                                                 *
43  ! nspeciesdim     = nspec for forward runs, 1 for backward runs              *
44  !                                                                            *
45  !*****************************************************************************
46
47  use unc_mod
48  use outg_mod
49  use par_mod
50  use com_mod
51
52  implicit none
53
54  integer :: itime,itage,i,ix,jy,ixp,jyp,kz,ks,n,nage
55  integer :: il,ind,indz,indzp,nrelpointer
56  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
57  real :: weight,hx,hy,hz,h,xd,yd,zd,xkern,r2,c(maxspec),ddx,ddy
58  real :: rhoprof(2),rhoi
59  real :: xl,yl,wx,wy,w
60  real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
61
62
63  ! For forward simulations, make a loop over the number of species;
64  ! for backward simulations, make an additional loop over the
65  ! releasepoints
66  !***************************************************************************
67
68
69  do i=1,numpart
70    if (itra1(i).ne.itime) goto 20
71
72  ! Determine age class of the particle
73    itage=abs(itra1(i)-itramem(i))
74    do nage=1,nageclass
75      if (itage.lt.lage(nage)) goto 33
76    end do
7733   continue
78
79
80  ! For special runs, interpolate the air density to the particle position
81  !************************************************************************
82  !***********************************************************************
83  !AF IND_SOURCE switches between different units for concentrations at the source
84  !Af    NOTE that in backward simulations the release of particles takes place
85  !Af    at the receptor and the sampling at the source.
86  !Af          1="mass"
87  !Af          2="mass mixing ratio"
88  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
89  !Af          1="mass"
90  !Af          2="mass mixing ratio"
91
92  !Af switches for the conccalcfile:
93  !AF IND_SAMP =  0 : xmass * 1
94  !Af IND_SAMP = -1 : xmass / rho
95
96  !Af ind_samp is defined in readcommand.f
97
98    if ( ind_samp .eq. -1 ) then
99
100      ix=int(xtra1(i))
101      jy=int(ytra1(i))
102      ixp=ix+1
103      jyp=jy+1
104      ddx=xtra1(i)-real(ix)
105      ddy=ytra1(i)-real(jy)
106      rddx=1.-ddx
107      rddy=1.-ddy
108      p1=rddx*rddy
109      p2=ddx*rddy
110      p3=rddx*ddy
111      p4=ddx*ddy
112
113      do il=2,nz
114        if (height(il).gt.ztra1(i)) then
115          indz=il-1
116          indzp=il
117          goto 6
118        endif
119      end do
1206     continue
121
122      dz1=ztra1(i)-height(indz)
123      dz2=height(indzp)-ztra1(i)
124      dz=1./(dz1+dz2)
125
126  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
127  !*****************************************************************************
128      do ind=indz,indzp
129        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
130             +p2*rho(ixp,jy ,ind,2) &
131             +p3*rho(ix ,jyp,ind,2) &
132             +p4*rho(ixp,jyp,ind,2)
133      end do
134      rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
135   elseif (ind_samp.eq.0) then
136      rhoi = 1.
137   endif
138
139
140  !****************************************************************************
141  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
142  !****************************************************************************
143
144
145  ! For backward simulations, look from which release point the particle comes from
146  ! For domain-filling trajectory option, npoint contains a consecutive particle
147  ! number, not the release point information. Therefore, nrelpointer is set to 1
148  ! for the domain-filling option.
149  !*****************************************************************************
150
151    if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
152       nrelpointer=1
153    else
154       nrelpointer=npoint(i)
155    endif
156
157    do kz=1,numzgrid                ! determine height of cell
158      if (outheight(kz).gt.ztra1(i)) goto 21
159    end do
16021   continue
161    if (kz.le.numzgrid) then           ! inside output domain
162
163
164  !********************************
165  ! Do everything for mother domain
166  !********************************
167
168      xl=(xtra1(i)*dx+xoutshift)/dxout
169      yl=(ytra1(i)*dy+youtshift)/dyout
170      ix=int(xl)
171      if (xl.lt.0.) ix=ix-1
172      jy=int(yl)
173      if (yl.lt.0.) jy=jy-1
174
175  ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
176
177
178  ! For particles aged less than 3 hours, attribute particle mass to grid cell
179  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
180  ! For older particles, use the uniform kernel.
181  ! If a particle is close to the domain boundary, do not use the kernel either.
182  !*****************************************************************************
183
184      if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
185           (xl.gt.real(numxgrid-1)-0.5).or. &
186           (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
187        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
188             (jy.le.numygrid-1)) then
189          do ks=1,nspec
190            gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
191                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
192                 xmass1(i,ks)/rhoi*weight
193          end do
194        endif
195
196      else                                 ! attribution via uniform kernel
197
198        ddx=xl-real(ix)                   ! distance to left cell border
199        ddy=yl-real(jy)                   ! distance to lower cell border
200        if (ddx.gt.0.5) then
201          ixp=ix+1
202          wx=1.5-ddx
203        else
204          ixp=ix-1
205          wx=0.5+ddx
206        endif
207
208        if (ddy.gt.0.5) then
209          jyp=jy+1
210          wy=1.5-ddy
211        else
212          jyp=jy-1
213          wy=0.5+ddy
214        endif
215
216
217  ! Determine mass fractions for four grid points
218  !**********************************************
219
220        if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
221          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
222            w=wx*wy
223            do ks=1,nspec
224              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
225                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
226                   xmass1(i,ks)/rhoi*weight*w
227            end do
228          endif
229
230          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
231            w=wx*(1.-wy)
232            do ks=1,nspec
233              gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
234                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
235                   xmass1(i,ks)/rhoi*weight*w
236            end do
237          endif
238        endif
239
240
241        if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
242          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
243            w=(1.-wx)*(1.-wy)
244            do ks=1,nspec
245              gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
246                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
247                   xmass1(i,ks)/rhoi*weight*w
248            end do
249          endif
250
251          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
252            w=(1.-wx)*wy
253            do ks=1,nspec
254              gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
255                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
256                   xmass1(i,ks)/rhoi*weight*w
257            end do
258          endif
259        endif
260      endif
261
262
263
264  !************************************
265  ! Do everything for the nested domain
266  !************************************
267
268      if (nested_output.eq.1) then
269        xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
270        yl=(ytra1(i)*dy+youtshiftn)/dyoutn
271        ix=int(xl)
272        if (xl.lt.0.) ix=ix-1
273        jy=int(yl)
274        if (yl.lt.0.) jy=jy-1
275
276
277  ! For particles aged less than 3 hours, attribute particle mass to grid cell
278  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
279  ! For older particles, use the uniform kernel.
280  ! If a particle is close to the domain boundary, do not use the kernel either.
281  !*****************************************************************************
282
283        if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
284             (xl.gt.real(numxgridn-1)-0.5).or. &
285             (yl.gt.real(numygridn-1)-0.5)) then             ! no kernel, direct attribution to grid cell
286          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
287               (jy.le.numygridn-1)) then
288            do ks=1,nspec
289              griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
290                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
291                   xmass1(i,ks)/rhoi*weight
292            end do
293          endif
294
295        else                                 ! attribution via uniform kernel
296
297          ddx=xl-real(ix)                   ! distance to left cell border
298          ddy=yl-real(jy)                   ! distance to lower cell border
299          if (ddx.gt.0.5) then
300            ixp=ix+1
301            wx=1.5-ddx
302          else
303            ixp=ix-1
304            wx=0.5+ddx
305          endif
306
307          if (ddy.gt.0.5) then
308            jyp=jy+1
309            wy=1.5-ddy
310          else
311            jyp=jy-1
312            wy=0.5+ddy
313          endif
314
315
316  ! Determine mass fractions for four grid points
317  !**********************************************
318
319          if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
320            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
321              w=wx*wy
322              do ks=1,nspec
323                griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
324                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
325                     xmass1(i,ks)/rhoi*weight*w
326              end do
327            endif
328
329            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
330              w=wx*(1.-wy)
331              do ks=1,nspec
332                griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
333                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
334                     xmass1(i,ks)/rhoi*weight*w
335              end do
336            endif
337          endif
338
339
340          if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
341            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
342              w=(1.-wx)*(1.-wy)
343              do ks=1,nspec
344                griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
345                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
346                     xmass1(i,ks)/rhoi*weight*w
347              end do
348            endif
349
350            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
351              w=(1.-wx)*wy
352              do ks=1,nspec
353                griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
354                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
355                     xmass1(i,ks)/rhoi*weight*w
356              end do
357            endif
358          endif
359        endif
360
361      endif
362    endif
36320  continue
364  end do
365
366  !***********************************************************************
367  ! 2. Evaluate concentrations at receptor points, using the kernel method
368  !***********************************************************************
369
370  do n=1,numreceptor
371
372
373  ! Reset concentrations
374  !*********************
375
376    do ks=1,nspec
377      c(ks)=0.
378    end do
379
380
381  ! Estimate concentration at receptor
382  !***********************************
383
384    do i=1,numpart
385
386      if (itra1(i).ne.itime) goto 40
387      itage=abs(itra1(i)-itramem(i))
388
389      hz=min(50.+0.3*sqrt(real(itage)),hzmax)
390      zd=ztra1(i)/hz
391      if (zd.gt.1.) goto 40          ! save computing time, leave loop
392
393      hx=min((0.29+2.222e-3*sqrt(real(itage)))*dx+ &
394           real(itage)*1.2e-5,hxmax)                     ! 80 km/day
395      xd=(xtra1(i)-xreceptor(n))/hx
396      if (xd*xd.gt.1.) goto 40       ! save computing time, leave loop
397
398      hy=min((0.18+1.389e-3*sqrt(real(itage)))*dy+ &
399           real(itage)*7.5e-6,hymax)                     ! 80 km/day
400      yd=(ytra1(i)-yreceptor(n))/hy
401      if (yd*yd.gt.1.) goto 40       ! save computing time, leave loop
402      h=hx*hy*hz
403
404      r2=xd*xd+yd*yd+zd*zd
405      if (r2.lt.1.) then
406        xkern=factor*(1.-r2)
407        do ks=1,nspec
408          c(ks)=c(ks)+xmass1(i,ks)*xkern/h
409        end do
410      endif
41140    continue
412    end do
413
414    do ks=1,nspec
415      creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n)
416    end do
417  end do
418
419end subroutine conccalc
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG