source: flexpart.git/src/conccalc.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was fe32dca, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Renamed lnokernel, corrected default setting.

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