conccalc_mpi.f90 18.1 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
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
43
44
  integer :: mind2
  ! mind2        eso: pointer to 2nd windfield in memory
45
46
47
48
49
50
51
52

  ! 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
53
54
  mind2=memind(2)

55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
  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

99
100
101
102
103
104
! eso: Temporary fix for particle exactly at north pole
      if (jyp >= nymax) then
        ! write(*,*) 'WARNING: conccalc.f90 jyp >= nymax'
        jyp=jyp-1
      end if

105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
      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
121
122
123
124
125
126
127
128
129
        ! 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)

130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
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
      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.
  !*****************************************************************************

181
182
      if ((.not.lusekerneloutput).or.(itage.lt.10800).or. &
           (xl.lt.0.5).or.(yl.lt.0.5).or. &
183
184
185
186
           (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
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
          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
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
        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
236
237
238
239
240
241
242
243
244
            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)= &
245
246
                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
247
248
               end do
            endif
249
250
251
252
          endif

          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
            w=wx*(1.-wy)
253
254
255
256
257
258
259
260
261
            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)= &
262
263
                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
264
265
               end do
             endif
266
267
268
269
270
271
272
          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)
273
274
275
276
277
278
279
280
281
            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)= &
282
283
                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
284
285
               end do
            endif
286
287
288
289
          endif

          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
            w=(1.-wx)*wy
290
291
292
293
294
295
296
297
298
            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)= &
299
300
                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
301
302
               end do
            endif
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
          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. &
330
331
             (yl.gt.real(numygridn-1)-0.5).or.((.not.lusekerneloutput))) then
! no kernel, direct attribution to grid cell
332
333
          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
               (jy.le.numygridn-1)) then
334
335
336
            if (DRYBKDEP.or.WETBKDEP) then
               do ks=1,nspec
                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
337
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
                   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
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
          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
383
384
385
386
387
388
389
390
391
              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)= &
392
393
                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
394
395
                 end do
              endif
396
397
398
399
            endif

            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
              w=wx*(1.-wy)
400
401
402
403
404
405
406
407
408
              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)= &
409
410
                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
411
412
                 end do
              endif
413
414
415
416
417
418
419
            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)
420
421
422
423
424
425
426
427
428
              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)= &
429
430
                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
431
432
                 end do
              endif
433
434
435
436
            endif

            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
              w=(1.-wx)*wy
437
438
439
440
441
442
443
444
445
              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)= &
446
447
                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
448
449
                 end do
              endif
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
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
            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