source: flexpart.git/src/conccalc.f90 @ 6985a98

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 6985a98 was 6985a98, checked in by Sabine <sabine.eckhardt@…>, 7 years ago

compiles after merge scavenging into test dev

  • Property mode set to 100644
File size: 18.4 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!  integer xscav_count
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!  xscav_count=0
68  do i=1,numpart
69    if (itra1(i).ne.itime) goto 20
70
71  ! Determine age class of the particle
72    itage=abs(itra1(i)-itramem(i))
73    do nage=1,nageclass
74      if (itage.lt.lage(nage)) goto 33
75    end do
7633   continue
77
78!  if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
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,memind(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
176
177  ! For particles aged less than 3 hours, attribute particle mass to grid cell
178  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
179  ! For older particles, use the uniform kernel.
180  ! If a particle is close to the domain boundary, do not use the kernel either.
181  !*****************************************************************************
182
183      if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
184           (xl.gt.real(numxgrid-1)-0.5).or. &
185           (yl.gt.real(numygrid-1)-0.5))) then             ! no kernel, direct attribution to grid cell
186        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
187             (jy.le.numygrid-1)) then
188          if (DRYBKDEP.or.WETBKDEP) 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*max(xscav_frac1(i,ks),0.0)
193             end do
194          else
195             do ks=1,nspec
196               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
197                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
198                 xmass1(i,ks)/rhoi*weight
199             end do
200          endif
201        endif
202
203      else                                 ! attribution via uniform kernel
204
205        ddx=xl-real(ix)                   ! distance to left cell border
206        ddy=yl-real(jy)                   ! distance to lower cell border
207        if (ddx.gt.0.5) then
208          ixp=ix+1
209          wx=1.5-ddx
210        else
211          ixp=ix-1
212          wx=0.5+ddx
213        endif
214
215        if (ddy.gt.0.5) then
216          jyp=jy+1
217          wy=1.5-ddy
218        else
219          jyp=jy-1
220          wy=0.5+ddy
221        endif
222
223
224  ! Determine mass fractions for four grid points
225  !**********************************************
226
227        if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
228          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
229            w=wx*wy
230            if (DRYBKDEP.or.WETBKDEP) then
231               do ks=1,nspec
232                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
233                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
234                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
235               end do
236            else
237               do ks=1,nspec
238                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
239                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
240                   xmass1(i,ks)/rhoi*weight*w
241               end do
242            endif
243          endif
244
245          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
246            w=wx*(1.-wy)
247            if (DRYBKDEP.or.WETBKDEP) then
248              do ks=1,nspec
249                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
250                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
251                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
252               end do
253             else
254              do ks=1,nspec
255                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
256                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
257                   xmass1(i,ks)/rhoi*weight*w
258               end do
259             endif
260          endif
261        endif !ix ge 0
262
263
264        if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
265          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
266            w=(1.-wx)*(1.-wy)
267            if (DRYBKDEP.or.WETBKDEP) then
268               do ks=1,nspec
269                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
270                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
271                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
272               end do
273            else
274               do ks=1,nspec
275                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
276                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
277                   xmass1(i,ks)/rhoi*weight*w
278               end do
279            endif
280          endif
281
282          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
283            w=(1.-wx)*wy
284            if (DRYBKDEP.or.WETBKDEP) then
285               do ks=1,nspec
286                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
287                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
288                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
289               end do
290            else
291               do ks=1,nspec
292                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
293                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
294                   xmass1(i,ks)/rhoi*weight*w
295               end do
296            endif
297          endif
298        endif !ixp ge 0
299     endif
300
301  !************************************
302  ! Do everything for the nested domain
303  !************************************
304
305      if (nested_output.eq.1) then
306        xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
307        yl=(ytra1(i)*dy+youtshiftn)/dyoutn
308        ix=int(xl)
309        if (xl.lt.0.) ix=ix-1
310        jy=int(yl)
311        if (yl.lt.0.) jy=jy-1
312
313
314  ! For particles aged less than 3 hours, attribute particle mass to grid cell
315  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
316  ! For older particles, use the uniform kernel.
317  ! If a particle is close to the domain boundary, do not use the kernel either.
318  !*****************************************************************************
319
320        if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
321             (xl.gt.real(numxgridn-1)-0.5).or. &
322             (yl.gt.real(numygridn-1)-0.5).or.(.not.usekernel)) then             ! no kernel, direct attribution to grid cell
323          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
324               (jy.le.numygridn-1)) then
325            if (DRYBKDEP.or.WETBKDEP) then
326               do ks=1,nspec
327                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
328                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
329                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
330               end do
331            else
332               do ks=1,nspec
333                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
334                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
335                   xmass1(i,ks)/rhoi*weight
336               end do
337            endif
338          endif
339
340        else                                 ! attribution via uniform kernel
341
342          ddx=xl-real(ix)                   ! distance to left cell border
343          ddy=yl-real(jy)                   ! distance to lower cell border
344          if (ddx.gt.0.5) then
345            ixp=ix+1
346            wx=1.5-ddx
347          else
348            ixp=ix-1
349            wx=0.5+ddx
350          endif
351
352          if (ddy.gt.0.5) then
353            jyp=jy+1
354            wy=1.5-ddy
355          else
356            jyp=jy-1
357            wy=0.5+ddy
358          endif
359
360
361  ! Determine mass fractions for four grid points
362  !**********************************************
363
364          if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
365            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
366              w=wx*wy
367              if (DRYBKDEP.or.WETBKDEP) then
368                 do ks=1,nspec
369                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
370                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
371                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
372                 end do
373              else
374                do ks=1,nspec
375                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
376                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
377                     xmass1(i,ks)/rhoi*weight*w
378                 end do
379              endif
380            endif
381
382            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
383              w=wx*(1.-wy)
384              if (DRYBKDEP.or.WETBKDEP) then
385                 do ks=1,nspec
386                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
387                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
388                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
389                 end do
390              else
391                 do ks=1,nspec
392                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
393                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
394                     xmass1(i,ks)/rhoi*weight*w
395                 end do
396              endif
397            endif
398          endif
399
400
401          if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
402            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
403              w=(1.-wx)*(1.-wy)
404              if (DRYBKDEP.or.WETBKDEP) then
405                 do ks=1,nspec
406                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
407                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
408                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
409                 end do
410              else
411                 do ks=1,nspec
412                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
413                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
414                     xmass1(i,ks)/rhoi*weight*w
415                 end do
416              endif
417            endif
418
419            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
420              w=(1.-wx)*wy
421              if (DRYBKDEP.or.WETBKDEP) then
422                 do ks=1,nspec
423                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
424                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
425                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
426                 end do
427              else
428                 do ks=1,nspec
429                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
430                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
431                     xmass1(i,ks)/rhoi*weight*w
432                 end do
433              endif
434            endif
435          endif
436        endif
437      endif
438    endif
43920  continue
440  end do
441!  write(*,*) 'xscav count:',xscav_count
442
443  !***********************************************************************
444  ! 2. Evaluate concentrations at receptor points, using the kernel method
445  !***********************************************************************
446
447  do n=1,numreceptor
448
449
450  ! Reset concentrations
451  !*********************
452
453    do ks=1,nspec
454      c(ks)=0.
455    end do
456
457
458  ! Estimate concentration at receptor
459  !***********************************
460
461    do i=1,numpart
462
463      if (itra1(i).ne.itime) goto 40
464      itage=abs(itra1(i)-itramem(i))
465
466      hz=min(50.+0.3*sqrt(real(itage)),hzmax)
467      zd=ztra1(i)/hz
468      if (zd.gt.1.) goto 40          ! save computing time, leave loop
469
470      hx=min((0.29+2.222e-3*sqrt(real(itage)))*dx+ &
471           real(itage)*1.2e-5,hxmax)                     ! 80 km/day
472      xd=(xtra1(i)-xreceptor(n))/hx
473      if (xd*xd.gt.1.) goto 40       ! save computing time, leave loop
474
475      hy=min((0.18+1.389e-3*sqrt(real(itage)))*dy+ &
476           real(itage)*7.5e-6,hymax)                     ! 80 km/day
477      yd=(ytra1(i)-yreceptor(n))/hy
478      if (yd*yd.gt.1.) goto 40       ! save computing time, leave loop
479      h=hx*hy*hz
480
481      r2=xd*xd+yd*yd+zd*zd
482      if (r2.lt.1.) then
483        xkern=factor*(1.-r2)
484        do ks=1,nspec
485          c(ks)=c(ks)+xmass1(i,ks)*xkern/h
486        end do
487      endif
48840    continue
489    end do
490
491    do ks=1,nspec
492      creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n)
493    end do
494  end do
495
496end subroutine conccalc
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG