conccalc_mpi.f90 19.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine conccalc(itime,weight)
  !                 i     i
  !*****************************************************************************
  !                                                                            *
  !     Calculation of the concentrations on a regular grid using volume       *
  !     sampling                                                               *
  !                                                                            *
  !     Author: A. Stohl                                                       *
  !                                                                            *
  !     24 May 1996                                                            *
  !                                                                            *
  !     April 2000: Update to calculate age spectra                            *
  !                 Bug fix to avoid negative conc. at the domain boundaries,  *
  !                 as suggested by Petra Seibert                              *
  !                                                                            *
  !     2 July 2002: re-order if-statements in order to optimize CPU time      *
  !                                                                            *
  !                                                                            *
  !*****************************************************************************
  !                                                                            *
  ! Variables:                                                                 *
  ! nspeciesdim     = nspec for forward runs, 1 for backward runs              *
  !                                                                            *
  !*****************************************************************************

  use unc_mod
  use outg_mod
  use par_mod
  use com_mod
  use mpi_mod

  implicit none

  integer,intent(in) :: itime
  integer :: itage,i,ix,jy,ixp,jyp,kz,ks,n,nage
  integer :: il,ind,indz,indzp,nrelpointer
  real,intent(in) :: weight
  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
  real :: hx,hy,hz,h,xd,yd,zd,xkern,r2,c(maxspec),ddx,ddy
  real :: rhoprof(2),rhoi
  real :: xl,yl,wx,wy,w
  real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
Espen Sollum's avatar
Espen Sollum committed
64
65
  integer :: mind2
  ! mind2        eso: pointer to 2nd windfield in memory
66
67
68
69
70
71
72
73

  ! For forward simulations, make a loop over the number of species;
  ! for backward simulations, make an additional loop over the
  ! releasepoints
  !***************************************************************************

  if (mp_measure_time) call mpif_mtime('conccalc',0)

Espen Sollum's avatar
Espen Sollum committed
74
75
  mind2=memind(2)

76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
  do i=1,numpart
    if (itra1(i).ne.itime) goto 20

  ! Determine age class of the particle
    itage=abs(itra1(i)-itramem(i))
    do nage=1,nageclass
      if (itage.lt.lage(nage)) goto 33
    end do
33   continue


  ! For special runs, interpolate the air density to the particle position
  !************************************************************************
  !***********************************************************************
  !AF IND_SOURCE switches between different units for concentrations at the source
  !Af    NOTE that in backward simulations the release of particles takes place
  !Af    at the receptor and the sampling at the source.
  !Af          1="mass"
  !Af          2="mass mixing ratio"
  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
  !Af          1="mass"
  !Af          2="mass mixing ratio"

  !Af switches for the conccalcfile:
  !AF IND_SAMP =  0 : xmass * 1
  !Af IND_SAMP = -1 : xmass / rho

  !Af ind_samp is defined in readcommand.f

    if ( ind_samp .eq. -1 ) then

      ix=int(xtra1(i))
      jy=int(ytra1(i))
      ixp=ix+1
      jyp=jy+1
      ddx=xtra1(i)-real(ix)
      ddy=ytra1(i)-real(jy)
      rddx=1.-ddx
      rddy=1.-ddy
      p1=rddx*rddy
      p2=ddx*rddy
      p3=rddx*ddy
      p4=ddx*ddy

120
121
122
123
124
125
! eso: Temporary fix for particle exactly at north pole
      if (jyp >= nymax) then
        ! write(*,*) 'WARNING: conccalc.f90 jyp >= nymax'
        jyp=jyp-1
      end if

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
      do il=2,nz
        if (height(il).gt.ztra1(i)) then
          indz=il-1
          indzp=il
          goto 6
        endif
      end do
6     continue

      dz1=ztra1(i)-height(indz)
      dz2=height(indzp)-ztra1(i)
      dz=1./(dz1+dz2)

  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
  !*****************************************************************************
      do ind=indz,indzp
Espen Sollum's avatar
Espen Sollum committed
142
143
144
145
146
147
148
149
150
        ! rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
        !      +p2*rho(ixp,jy ,ind,2) &
        !      +p3*rho(ix ,jyp,ind,2) &
        !      +p4*rho(ixp,jyp,ind,2)
        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,mind2) &
             +p2*rho(ixp,jy ,ind,mind2) &
             +p3*rho(ix ,jyp,ind,mind2) &
             +p4*rho(ixp,jyp,ind,mind2)

151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
      end do
      rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
   elseif (ind_samp.eq.0) then
      rhoi = 1.
   endif


  !****************************************************************************
  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
  !****************************************************************************


  ! For backward simulations, look from which release point the particle comes from
  ! For domain-filling trajectory option, npoint contains a consecutive particle
  ! number, not the release point information. Therefore, nrelpointer is set to 1
  ! for the domain-filling option.
  !*****************************************************************************

    if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
       nrelpointer=1
    else
       nrelpointer=npoint(i)
    endif

    do kz=1,numzgrid                ! determine height of cell
      if (outheight(kz).gt.ztra1(i)) goto 21
    end do
21   continue
    if (kz.le.numzgrid) then           ! inside output domain


  !********************************
  ! Do everything for mother domain
  !********************************

      xl=(xtra1(i)*dx+xoutshift)/dxout
      yl=(ytra1(i)*dy+youtshift)/dyout
      ix=int(xl)
      if (xl.lt.0.) ix=ix-1
      jy=int(yl)
      if (yl.lt.0.) jy=jy-1

  ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl


  ! For particles aged less than 3 hours, attribute particle mass to grid cell
  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
  ! For older particles, use the uniform kernel.
  ! If a particle is close to the domain boundary, do not use the kernel either.
  !*****************************************************************************

202
      if (lnokernel.or.(itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
203
204
205
206
           (xl.gt.real(numxgrid-1)-0.5).or. &
           (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
             (jy.le.numygrid-1)) then
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
          if (DRYBKDEP.or.WETBKDEP) then
            do ks=1,nspec
              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
            end do
          else
            if (lparticlecountoutput) then
              do ks=1,nspec
                gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                     gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
              end do
            else
              do ks=1,nspec
                gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                     gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight
              end do
            end if
          endif
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
        endif

      else                                 ! attribution via uniform kernel

        ddx=xl-real(ix)                   ! distance to left cell border
        ddy=yl-real(jy)                   ! distance to lower cell border
        if (ddx.gt.0.5) then
          ixp=ix+1
          wx=1.5-ddx
        else
          ixp=ix-1
          wx=0.5+ddx
        endif

        if (ddy.gt.0.5) then
          jyp=jy+1
          wy=1.5-ddy
        else
          jyp=jy-1
          wy=0.5+ddy
        endif


  ! Determine mass fractions for four grid points
  !**********************************************

        if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
            w=wx*wy
256
257
258
259
260
261
262
263
264
            if (DRYBKDEP.or.WETBKDEP) then
               do ks=1,nspec
                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
               end do
            else
               do ks=1,nspec
                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
265
266
                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
267
268
               end do
            endif
269
270
271
272
          endif

          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
            w=wx*(1.-wy)
273
274
275
276
277
278
279
280
281
            if (DRYBKDEP.or.WETBKDEP) then
              do ks=1,nspec
                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
               end do
             else
              do ks=1,nspec
                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
282
283
                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
284
285
               end do
             endif
286
287
288
289
290
291
292
          endif
        endif


        if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
            w=(1.-wx)*(1.-wy)
293
294
295
296
297
298
299
300
301
            if (DRYBKDEP.or.WETBKDEP) then
               do ks=1,nspec
                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
               end do
            else
               do ks=1,nspec
                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
302
303
                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
304
305
               end do
            endif
306
307
308
309
          endif

          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
            w=(1.-wx)*wy
310
311
312
313
314
315
316
317
318
            if (DRYBKDEP.or.WETBKDEP) then
               do ks=1,nspec
                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
               end do
            else
               do ks=1,nspec
                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
319
320
                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
321
322
               end do
            endif
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
          endif
        endif
      endif



  !************************************
  ! Do everything for the nested domain
  !************************************

      if (nested_output.eq.1) then
        xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
        yl=(ytra1(i)*dy+youtshiftn)/dyoutn
        ix=int(xl)
        if (xl.lt.0.) ix=ix-1
        jy=int(yl)
        if (yl.lt.0.) jy=jy-1


  ! For particles aged less than 3 hours, attribute particle mass to grid cell
  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
  ! For older particles, use the uniform kernel.
  ! If a particle is close to the domain boundary, do not use the kernel either.
  !*****************************************************************************

        if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
             (xl.gt.real(numxgridn-1)-0.5).or. &
             (yl.gt.real(numygridn-1)-0.5)) then             ! no kernel, direct attribution to grid cell
          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
               (jy.le.numygridn-1)) then
353
354
355
            if (DRYBKDEP.or.WETBKDEP) then
               do ks=1,nspec
                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
356
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
               end do
            else
              if (lparticlecountoutput) then
                do ks=1,nspec
                  griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                       griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+1
                end do
              else
                do ks=1,nspec
                  griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                       griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                       xmass1(i,ks)/rhoi*weight
                end do
              endif
            endif
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
          endif

        else                                 ! attribution via uniform kernel

          ddx=xl-real(ix)                   ! distance to left cell border
          ddy=yl-real(jy)                   ! distance to lower cell border
          if (ddx.gt.0.5) then
            ixp=ix+1
            wx=1.5-ddx
          else
            ixp=ix-1
            wx=0.5+ddx
          endif

          if (ddy.gt.0.5) then
            jyp=jy+1
            wy=1.5-ddy
          else
            jyp=jy-1
            wy=0.5+ddy
          endif


  ! Determine mass fractions for four grid points
  !**********************************************

          if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
              w=wx*wy
402
403
404
405
406
407
408
409
410
              if (DRYBKDEP.or.WETBKDEP) then
                 do ks=1,nspec
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
                 end do
              else
                do ks=1,nspec
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
411
412
                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
413
414
                 end do
              endif
415
416
417
418
            endif

            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
              w=wx*(1.-wy)
419
420
421
422
423
424
425
426
427
              if (DRYBKDEP.or.WETBKDEP) then
                 do ks=1,nspec
                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
                 end do
              else
                 do ks=1,nspec
                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
428
429
                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
430
431
                 end do
              endif
432
433
434
435
436
437
438
            endif
          endif


          if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
              w=(1.-wx)*(1.-wy)
439
440
441
442
443
444
445
446
447
              if (DRYBKDEP.or.WETBKDEP) then
                 do ks=1,nspec
                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
                 end do
              else
                 do ks=1,nspec
                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
448
449
                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
450
451
                 end do
              endif
452
453
454
455
            endif

            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
              w=(1.-wx)*wy
456
457
458
459
460
461
462
463
464
              if (DRYBKDEP.or.WETBKDEP) then
                 do ks=1,nspec
                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
                 end do
              else
                 do ks=1,nspec
                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
465
466
                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
467
468
                 end do
              endif
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
            endif
          endif
        endif

      endif
    endif
20  continue
  end do

  !***********************************************************************
  ! 2. Evaluate concentrations at receptor points, using the kernel method
  !***********************************************************************

  do n=1,numreceptor


  ! Reset concentrations
  !*********************

    do ks=1,nspec
      c(ks)=0.
    end do


  ! Estimate concentration at receptor
  !***********************************

    do i=1,numpart

      if (itra1(i).ne.itime) goto 40
      itage=abs(itra1(i)-itramem(i))

      hz=min(50.+0.3*sqrt(real(itage)),hzmax)
      zd=ztra1(i)/hz
      if (zd.gt.1.) goto 40          ! save computing time, leave loop

      hx=min((0.29+2.222e-3*sqrt(real(itage)))*dx+ &
           real(itage)*1.2e-5,hxmax)                     ! 80 km/day
      xd=(xtra1(i)-xreceptor(n))/hx
      if (xd*xd.gt.1.) goto 40       ! save computing time, leave loop

      hy=min((0.18+1.389e-3*sqrt(real(itage)))*dy+ &
           real(itage)*7.5e-6,hymax)                     ! 80 km/day
      yd=(ytra1(i)-yreceptor(n))/hy
      if (yd*yd.gt.1.) goto 40       ! save computing time, leave loop
      h=hx*hy*hz

      r2=xd*xd+yd*yd+zd*zd
      if (r2.lt.1.) then
        xkern=factor*(1.-r2)
        do ks=1,nspec
          c(ks)=c(ks)+xmass1(i,ks)*xkern/h
        end do
      endif
40    continue
    end do

    do ks=1,nspec
      creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n)
    end do
  end do

  if (mp_measure_time) call mpif_mtime('conccalc',1)

end subroutine conccalc