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

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

Renamed lnokernel, corrected default setting.

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