conccalc.f90 17.8 KB
Newer Older
Matthias Langer's avatar
 
Matthias Langer committed
1
subroutine conccalc(itime,weight)
2
  !                      i     i
Matthias Langer's avatar
 
Matthias Langer committed
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
  !*****************************************************************************
  !                                                                            *
  !     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

  implicit none

  integer :: itime,itage,i,ix,jy,ixp,jyp,kz,ks,n,nage
  integer :: il,ind,indz,indzp,nrelpointer
  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
  real :: weight,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.
Sabine's avatar
Sabine committed
40
!  integer xscav_count
Matthias Langer's avatar
 
Matthias Langer committed
41
42
43
44
45

  ! For forward simulations, make a loop over the number of species;
  ! for backward simulations, make an additional loop over the
  ! releasepoints
  !***************************************************************************
Sabine's avatar
Sabine committed
46
!  xscav_count=0
Matthias Langer's avatar
 
Matthias Langer committed
47
48
49
50
51
52
53
54
55
56
  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

Sabine's avatar
Sabine committed
57
!  if (xscav_frac1(i,1).lt.0) xscav_count=xscav_count+1
58
           
Matthias Langer's avatar
 
Matthias Langer committed
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
  ! 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

92
93
94
95
96
97
! eso: Temporary fix for particle exactly at north pole
      if (jyp >= nymax) then
      !  write(*,*) 'WARNING: conccalc.f90 jyp >= nymax'
        jyp=jyp-1
      end if

Matthias Langer's avatar
 
Matthias Langer committed
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
      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
114
        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,memind(2)) &
Matthias Langer's avatar
 
Matthias Langer committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
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
             +p2*rho(ixp,jy ,ind,2) &
             +p3*rho(ix ,jyp,ind,2) &
             +p4*rho(ixp,jyp,ind,2)
      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



  ! 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.
  !*****************************************************************************

168
169
      if ((.not.lusekerneloutput).or.(itage.lt.10800).or. &
           (xl.lt.0.5).or.(yl.lt.0.5).or. &
Matthias Langer's avatar
 
Matthias Langer committed
170
           (xl.gt.real(numxgrid-1)-0.5).or. &
171
           (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
Matthias Langer's avatar
 
Matthias Langer committed
172
173
        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
             (jy.le.numygrid-1)) then
174
          if (DRYBKDEP.or.WETBKDEP) then
175
176
177
178
179
            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
180
          else
181
182
183
184
185
186
187
188
189
190
191
192
            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
193
          endif
Matthias Langer's avatar
 
Matthias Langer committed
194
195
        endif

Sabine's avatar
Sabine committed
196
      else                                 ! attribution via uniform kernel 
Matthias Langer's avatar
 
Matthias Langer committed
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222

        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
223
            if (DRYBKDEP.or.WETBKDEP) then
224
225
226
               do ks=1,nspec
                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
227
                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
228
229
230
231
               end do
            else
               do ks=1,nspec
                 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
232
233
                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
234
235
               end do
            endif
Matthias Langer's avatar
 
Matthias Langer committed
236
237
238
239
          endif

          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
            w=wx*(1.-wy)
240
            if (DRYBKDEP.or.WETBKDEP) then
241
242
243
              do ks=1,nspec
                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
244
                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
245
246
247
248
               end do
             else
              do ks=1,nspec
                 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
249
250
                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
251
252
               end do
             endif
Matthias Langer's avatar
 
Matthias Langer committed
253
          endif
254
        endif !ix ge 0
Matthias Langer's avatar
 
Matthias Langer committed
255
256
257
258
259


        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)
260
            if (DRYBKDEP.or.WETBKDEP) then
261
262
263
               do ks=1,nspec
                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
264
                   xmass1(i,ks)/rhoi*w*weight*max(xscav_frac1(i,ks),0.0)
265
266
267
268
               end do
            else
               do ks=1,nspec
                 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
269
270
                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
271
272
               end do
            endif
Matthias Langer's avatar
 
Matthias Langer committed
273
274
275
276
          endif

          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
            w=(1.-wx)*wy
277
            if (DRYBKDEP.or.WETBKDEP) then
278
279
280
               do ks=1,nspec
                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
281
                   xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
282
283
284
285
               end do
            else
               do ks=1,nspec
                 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
286
287
                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                   xmass1(i,ks)/rhoi*weight*w
288
289
               end do
            endif
Matthias Langer's avatar
 
Matthias Langer committed
290
          endif
291
292
        endif !ixp ge 0
     endif
Matthias Langer's avatar
 
Matthias Langer committed
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314

  !************************************
  ! 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. &
315
316
             (yl.gt.real(numygridn-1)-0.5).or.((.not.lusekerneloutput))) then
! no kernel, direct attribution to grid cell
Matthias Langer's avatar
 
Matthias Langer committed
317
318
          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
               (jy.le.numygridn-1)) then
319
            if (DRYBKDEP.or.WETBKDEP) then
320
321
322
               do ks=1,nspec
                 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
323
                   xmass1(i,ks)/rhoi*weight*max(xscav_frac1(i,ks),0.0)
324
325
               end do
            else
326
327
328
329
330
331
332
333
334
335
336
337
              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
338
            endif
Matthias Langer's avatar
 
Matthias Langer committed
339
          endif
340
          
Matthias Langer's avatar
 
Matthias Langer committed
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
        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
368
              if (DRYBKDEP.or.WETBKDEP) then
369
370
371
                 do ks=1,nspec
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
372
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
373
374
375
376
                 end do
              else
                do ks=1,nspec
                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
377
378
                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
379
380
                 end do
              endif
Matthias Langer's avatar
 
Matthias Langer committed
381
382
383
384
            endif

            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
              w=wx*(1.-wy)
385
              if (DRYBKDEP.or.WETBKDEP) then
386
387
388
                 do ks=1,nspec
                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
389
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
390
391
392
393
                 end do
              else
                 do ks=1,nspec
                   griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
394
395
                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
396
397
                 end do
              endif
Matthias Langer's avatar
 
Matthias Langer committed
398
399
400
401
402
403
404
            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)
405
              if (DRYBKDEP.or.WETBKDEP) then
406
407
408
                 do ks=1,nspec
                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
409
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
410
411
412
413
                 end do
              else
                 do ks=1,nspec
                   griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
414
415
                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
416
417
                 end do
              endif
Matthias Langer's avatar
 
Matthias Langer committed
418
419
420
421
            endif

            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
              w=(1.-wx)*wy
422
              if (DRYBKDEP.or.WETBKDEP) then
423
424
425
                 do ks=1,nspec
                   griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
426
                     xmass1(i,ks)/rhoi*weight*w*max(xscav_frac1(i,ks),0.0)
427
428
429
430
                 end do
              else
                 do ks=1,nspec
                    griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
Matthias Langer's avatar
 
Matthias Langer committed
431
432
                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
                     xmass1(i,ks)/rhoi*weight*w
433
434
                 end do
              endif
Matthias Langer's avatar
 
Matthias Langer committed
435
436
437
438
439
440
441
            endif
          endif
        endif
      endif
    endif
20  continue
  end do
Sabine's avatar
Sabine committed
442
!  write(*,*) 'xscav count:',xscav_count
Matthias Langer's avatar
 
Matthias Langer committed
443
444
445
446
447
448
449
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

  !***********************************************************************
  ! 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

end subroutine conccalc