source: flexpart.git/src/conccalc.f90 @ 278f4ed

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

message on kernel use in timemanager

  • Property mode set to 100644
File size: 18.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 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  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
79  ! For special runs, interpolate the air density to the particle position
80  !************************************************************************
81  !***********************************************************************
82  !AF IND_SOURCE switches between different units for concentrations at the source
83  !Af    NOTE that in backward simulations the release of particles takes place
84  !Af    at the receptor and the sampling at the source.
85  !Af          1="mass"
86  !Af          2="mass mixing ratio"
87  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
88  !Af          1="mass"
89  !Af          2="mass mixing ratio"
90
91  !Af switches for the conccalcfile:
92  !AF IND_SAMP =  0 : xmass * 1
93  !Af IND_SAMP = -1 : xmass / rho
94
95  !Af ind_samp is defined in readcommand.f
96
97    if ( ind_samp .eq. -1 ) then
98
99      ix=int(xtra1(i))
100      jy=int(ytra1(i))
101      ixp=ix+1
102      jyp=jy+1
103      ddx=xtra1(i)-real(ix)
104      ddy=ytra1(i)-real(jy)
105      rddx=1.-ddx
106      rddy=1.-ddy
107      p1=rddx*rddy
108      p2=ddx*rddy
109      p3=rddx*ddy
110      p4=ddx*ddy
111
112      do il=2,nz
113        if (height(il).gt.ztra1(i)) then
114          indz=il-1
115          indzp=il
116          goto 6
117        endif
118      end do
1196     continue
120
121      dz1=ztra1(i)-height(indz)
122      dz2=height(indzp)-ztra1(i)
123      dz=1./(dz1+dz2)
124
125  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
126  !*****************************************************************************
127      do ind=indz,indzp
128        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
129             +p2*rho(ixp,jy ,ind,2) &
130             +p3*rho(ix ,jyp,ind,2) &
131             +p4*rho(ixp,jyp,ind,2)
132      end do
133      rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
134   elseif (ind_samp.eq.0) then
135      rhoi = 1.
136   endif
137
138
139  !****************************************************************************
140  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
141  !****************************************************************************
142
143
144  ! For backward simulations, look from which release point the particle comes from
145  ! For domain-filling trajectory option, npoint contains a consecutive particle
146  ! number, not the release point information. Therefore, nrelpointer is set to 1
147  ! for the domain-filling option.
148  !*****************************************************************************
149
150    if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
151       nrelpointer=1
152    else
153       nrelpointer=npoint(i)
154    endif
155
156    do kz=1,numzgrid                ! determine height of cell
157      if (outheight(kz).gt.ztra1(i)) goto 21
158    end do
15921   continue
160    if (kz.le.numzgrid) then           ! inside output domain
161
162
163  !********************************
164  ! Do everything for mother domain
165  !********************************
166
167      xl=(xtra1(i)*dx+xoutshift)/dxout
168      yl=(ytra1(i)*dy+youtshift)/dyout
169      ix=int(xl)
170      if (xl.lt.0.) ix=ix-1
171      jy=int(yl)
172      if (yl.lt.0.) jy=jy-1
173
174
175
176  ! For particles aged less than 3 hours, attribute particle mass to grid cell
177  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
178  ! For older particles, use the uniform kernel.
179  ! If a particle is close to the domain boundary, do not use the kernel either.
180  !*****************************************************************************
181
182      if (((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
183           (xl.gt.real(numxgrid-1)-0.5).or. &
184           (yl.gt.real(numygrid-1)-0.5)).or.(.not.usekernel)) then             ! no kernel, direct attribution to grid cell
185        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
186             (jy.le.numygrid-1)) then
187          if (DRYBKDEP.or.WETBKDEP) then
188             do ks=1,nspec
189               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
190                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
191                 xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
192             end do
193          else
194             do ks=1,nspec
195               gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
196                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
197                 xmass1(i,ks)/rhoi*weight
198             end do
199          endif
200        endif
201
202      else                                 ! attribution via uniform kernel
203
204        ddx=xl-real(ix)                   ! distance to left cell border
205        ddy=yl-real(jy)                   ! distance to lower cell border
206        if (ddx.gt.0.5) then
207          ixp=ix+1
208          wx=1.5-ddx
209        else
210          ixp=ix-1
211          wx=0.5+ddx
212        endif
213
214        if (ddy.gt.0.5) then
215          jyp=jy+1
216          wy=1.5-ddy
217        else
218          jyp=jy-1
219          wy=0.5+ddy
220        endif
221
222
223  ! Determine mass fractions for four grid points
224  !**********************************************
225
226        if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
227          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
228            w=wx*wy
229            if (DRYBKDEP.or.WETBKDEP) then
230               do ks=1,nspec
231                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
232                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
233                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
234               end do
235            else
236               do ks=1,nspec
237                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
238                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
239                   xmass1(i,ks)/rhoi*weight*w
240               end do
241            endif
242          endif
243
244          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
245            w=wx*(1.-wy)
246            if (DRYBKDEP.or.WETBKDEP) then
247              do ks=1,nspec
248                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
249                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
250                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
251               end do
252             else
253              do ks=1,nspec
254                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
255                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
256                   xmass1(i,ks)/rhoi*weight*w
257               end do
258             endif
259          endif
260        endif !ix ge 0
261
262
263        if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
264          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
265            w=(1.-wx)*(1.-wy)
266            if (DRYBKDEP.or.WETBKDEP) then
267               do ks=1,nspec
268                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
269                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
270                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
271               end do
272            else
273               do ks=1,nspec
274                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
275                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
276                   xmass1(i,ks)/rhoi*weight*w
277               end do
278            endif
279          endif
280
281          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
282            w=(1.-wx)*wy
283            if (DRYBKDEP.or.WETBKDEP) then
284               do ks=1,nspec
285                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
286                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
287                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
288               end do
289            else
290               do ks=1,nspec
291                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
292                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
293                   xmass1(i,ks)/rhoi*weight*w
294               end do
295            endif
296          endif
297        endif !ixp ge 0
298     endif
299
300  !************************************
301  ! Do everything for the nested domain
302  !************************************
303
304      if (nested_output.eq.1) then
305        xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
306        yl=(ytra1(i)*dy+youtshiftn)/dyoutn
307        ix=int(xl)
308        if (xl.lt.0.) ix=ix-1
309        jy=int(yl)
310        if (yl.lt.0.) jy=jy-1
311
312
313  ! For particles aged less than 3 hours, attribute particle mass to grid cell
314  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
315  ! For older particles, use the uniform kernel.
316  ! If a particle is close to the domain boundary, do not use the kernel either.
317  !*****************************************************************************
318
319        if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
320             (xl.gt.real(numxgridn-1)-0.5).or. &
321             (yl.gt.real(numygridn-1)-0.5).or.(.not.usekernel)) then             ! no kernel, direct attribution to grid cell
322          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
323               (jy.le.numygridn-1)) then
324            if (DRYBKDEP.or.WETBKDEP) then
325               do ks=1,nspec
326                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
327                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
328                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
329               end do
330            else
331               do ks=1,nspec
332                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
333                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
334                   xmass1(i,ks)/rhoi*weight
335               end do
336            endif
337          endif
338
339        else                                 ! attribution via uniform kernel
340
341          ddx=xl-real(ix)                   ! distance to left cell border
342          ddy=yl-real(jy)                   ! distance to lower cell border
343          if (ddx.gt.0.5) then
344            ixp=ix+1
345            wx=1.5-ddx
346          else
347            ixp=ix-1
348            wx=0.5+ddx
349          endif
350
351          if (ddy.gt.0.5) then
352            jyp=jy+1
353            wy=1.5-ddy
354          else
355            jyp=jy-1
356            wy=0.5+ddy
357          endif
358
359
360  ! Determine mass fractions for four grid points
361  !**********************************************
362
363          if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
364            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
365              w=wx*wy
366              if (DRYBKDEP.or.WETBKDEP) then
367                 do ks=1,nspec
368                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
369                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
370                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
371                 end do
372              else
373                do ks=1,nspec
374                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
375                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
376                     xmass1(i,ks)/rhoi*weight*w
377                 end do
378              endif
379            endif
380
381            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
382              w=wx*(1.-wy)
383              if (DRYBKDEP.or.WETBKDEP) then
384                 do ks=1,nspec
385                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
386                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
387                     xmass1(i,ks)/rhoi*weight*w*xscav_frac1(i,ks)
388                 end do
389              else
390                 do ks=1,nspec
391                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
392                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
393                     xmass1(i,ks)/rhoi*weight*w
394                 end do
395              endif
396            endif
397          endif
398
399
400          if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
401            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
402              w=(1.-wx)*(1.-wy)
403              if (DRYBKDEP.or.WETBKDEP) then
404                 do ks=1,nspec
405                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
406                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
407                     xmass1(i,ks)/rhoi*weight*w*xscav_frac1(i,ks)
408                 end do
409              else
410                 do ks=1,nspec
411                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
412                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
413                     xmass1(i,ks)/rhoi*weight*w
414                 end do
415              endif
416            endif
417
418            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
419              w=(1.-wx)*wy
420              if (DRYBKDEP.or.WETBKDEP) then
421                 do ks=1,nspec
422                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
423                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
424                     xmass1(i,ks)/rhoi*weight*w*xscav_frac1(i,ks)
425                 end do
426              else
427                 do ks=1,nspec
428                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
429                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
430                     xmass1(i,ks)/rhoi*weight*w
431                 end do
432              endif
433            endif
434          endif
435        endif
436      endif
437    endif
43820  continue
439  end do
440
441  !***********************************************************************
442  ! 2. Evaluate concentrations at receptor points, using the kernel method
443  !***********************************************************************
444
445  do n=1,numreceptor
446
447
448  ! Reset concentrations
449  !*********************
450
451    do ks=1,nspec
452      c(ks)=0.
453    end do
454
455
456  ! Estimate concentration at receptor
457  !***********************************
458
459    do i=1,numpart
460
461      if (itra1(i).ne.itime) goto 40
462      itage=abs(itra1(i)-itramem(i))
463
464      hz=min(50.+0.3*sqrt(real(itage)),hzmax)
465      zd=ztra1(i)/hz
466      if (zd.gt.1.) goto 40          ! save computing time, leave loop
467
468      hx=min((0.29+2.222e-3*sqrt(real(itage)))*dx+ &
469           real(itage)*1.2e-5,hxmax)                     ! 80 km/day
470      xd=(xtra1(i)-xreceptor(n))/hx
471      if (xd*xd.gt.1.) goto 40       ! save computing time, leave loop
472
473      hy=min((0.18+1.389e-3*sqrt(real(itage)))*dy+ &
474           real(itage)*7.5e-6,hymax)                     ! 80 km/day
475      yd=(ytra1(i)-yreceptor(n))/hy
476      if (yd*yd.gt.1.) goto 40       ! save computing time, leave loop
477      h=hx*hy*hz
478
479      r2=xd*xd+yd*yd+zd*zd
480      if (r2.lt.1.) then
481        xkern=factor*(1.-r2)
482        do ks=1,nspec
483          c(ks)=c(ks)+xmass1(i,ks)*xkern/h
484        end do
485      endif
48640    continue
487    end do
488
489    do ks=1,nspec
490      creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n)
491    end do
492  end do
493
494end subroutine conccalc
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG