source: branches/jerome/src_flexwrf_v3.1/conccalc_irreg.f90 @ 16

Last change on this file since 16 was 16, checked in by jebri, 11 years ago

sources for flexwrf v3.1

File size: 15.2 KB
Line 
1!***********************************************************************
2!* Copyright 2012,2013                                                *
3!* Jerome Brioude, Delia Arnold, Andreas Stohl, Wayne Angevine,       *
4!* John Burkhart, Massimo Cassiani, Adam Dingwell, Richard C Easter, Sabine Eckhardt,*
5!* Stephanie Evan, Jerome D Fast, Don Morton, Ignacio Pisso,          *
6!* Petra Seibert, Gerard Wotawa, Caroline Forster, Harald Sodemann,   *
7!*                                                                     *
8!* This file is part of FLEXPART WRF                                   *
9!*                                                                     *
10!* FLEXPART is free software: you can redistribute it and/or modify    *
11!* it under the terms of the GNU General Public License as published by*
12!* the Free Software Foundation, either version 3 of the License, or   *
13!* (at your option) any later version.                                 *
14!*                                                                     *
15!* FLEXPART is distributed in the hope that it will be useful,         *
16!* but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17!* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18!* GNU General Public License for more details.                        *
19!*                                                                     *
20!* You should have received a copy of the GNU General Public License   *
21!* along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
22!***********************************************************************
23
24subroutine conccalc_irreg(itime,weight)
25  !                      i     i
26  !*****************************************************************************
27  !                                                                            *
28  !     Calculation of the concentrations on a regular grid using volume       *
29  !     sampling                                                               *
30  !                                                                            *
31  !     Author: A. Stohl                                                       *
32  !                                                                            *
33  !     24 May 1996                                                            *
34  !                                                                            *
35  !     April 2000: Update to calculate age spectra                            *
36  !                 Bug fix to avoid negative conc. at the domain boundaries,  *
37  !                 as suggested by Petra Seibert                              *
38  !                                                                            *
39  !     2 July 2002: re-order if-statements in order to optimize CPU time      *
40  !                                                                            *
41  !                                                                            *
42  !*****************************************************************************
43  !                                                                            *
44  ! Variables:                                                                 *
45  ! nspeciesdim     = nspec for forward runs, 1 for backward runs              *
46  !                                                                            *
47  !*****************************************************************************
48
49  use unc_mod
50  use outg_mod
51  use par_mod
52  use com_mod
53
54  implicit none
55
56  integer :: itime,itage,i,ix,jy,ixp,jyp,kz,ks,n,nage
57  integer :: il,ind,indz,indzp,nrelpointer,ii,temp,incl
58  real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
59  real :: weight,hx,hy,hz,h,xd,yd,zd,xkern,r2,c(maxspec),ddx,ddy
60  real :: rhoprof(2),rhoi
61  real :: xl,yl,wx,wy,w
62  real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
63
64
65  ! For forward simulations, make a loop over the number of species;
66  ! for backward simulations, make an additional loop over the
67  ! releasepoints
68  !***************************************************************************
69
70
71  do i=1,numpart
72    if (itra1(i).ne.itime .or. abs(npoint(i)).gt.numpoint) goto 20
73
74  ! Determine age class of the particle
75    itage=abs(itra1(i)-itramem(i))
76    do nage=1,nageclass
77      if (itage.lt.lage(nage)) goto 33
78    end do
7933   continue
80
81
82  ! For special runs, interpolate the air density to the particle position
83  !************************************************************************
84  !***********************************************************************
85  !AF IND_SOURCE switches between different units for concentrations at the source
86  !Af    NOTE that in backward simulations the release of particles takes place
87  !Af    at the receptor and the sampling at the source.
88  !Af          1="mass"
89  !Af          2="mass mixing ratio"
90  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
91  !Af          1="mass"
92  !Af          2="mass mixing ratio"
93
94  !Af switches for the conccalcfile:
95  !AF IND_SAMP =  0 : xmass * 1
96  !Af IND_SAMP = -1 : xmass / rho
97
98  !Af ind_samp is defined in readcommand.f
99
100    if ( ind_samp .eq. -1 ) then
101
102      ix=int(xtra1(i))
103      jy=int(ytra1(i))
104      ixp=ix+1
105      jyp=jy+1
106      ddx=xtra1(i)-real(ix)
107      ddy=ytra1(i)-real(jy)
108      rddx=1.-ddx
109      rddy=1.-ddy
110      p1=rddx*rddy
111      p2=ddx*rddy
112      p3=rddx*ddy
113      p4=ddx*ddy
114
115      do il=2,nz
116        if (height(il).gt.ztra1(i)) then
117          indz=il-1
118          indzp=il
119          goto 6
120        endif
121      end do
1226     continue
123
124      dz1=ztra1(i)-height(indz)
125      dz2=height(indzp)-ztra1(i)
126      dz=1./(dz1+dz2)
127
128  ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed)
129  !*****************************************************************************
130       if (xtra1(i).ne.xtra1(i))  print*,i,itra1(i),xtra1(i),ytra1(i),ix,jy
131      do ind=indz,indzp
132        rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
133             +p2*rho(ixp,jy ,ind,2) &
134             +p3*rho(ix ,jyp,ind,2) &
135             +p4*rho(ixp,jyp,ind,2)
136      end do
137      rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
138   elseif (ind_samp.eq.0) then
139      rhoi = 1.
140   endif
141
142
143  !****************************************************************************
144  ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy
145  !****************************************************************************
146
147
148  ! For backward simulations, look from which release point the particle comes from
149  ! For domain-filling trajectory option, npoint contains a consecutive particle
150  ! number, not the release point information. Therefore, nrelpointer is set to 1
151  ! for the domain-filling option.
152  !*****************************************************************************
153
154    if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then
155       nrelpointer=1
156    else
157       nrelpointer=npoint(i)
158    endif
159
160    do kz=1,numzgrid                ! determine height of cell
161      if (outheight(kz).gt.ztra1(i)) goto 21
162    end do
16321   continue
164    if (kz.le.numzgrid) then           ! inside output domain
165
166
167  !********************************
168  ! Do everything for mother domain
169  !********************************
170
171      xl=(xtra1(i)*dx+xoutshift)/dxout
172      yl=(ytra1(i)*dy+youtshift)/dyout
173      ix=int(xl)
174      if (xl.lt.0.) ix=ix-1
175      jy=int(yl)
176      if (yl.lt.0.) jy=jy-1
177
178  ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl
179
180
181  ! For particles aged less than 3 hours, attribute particle mass to grid cell
182  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
183  ! For older particles, use the uniform kernel.
184  ! If a particle is close to the domain boundary, do not use the kernel either.
185  !*****************************************************************************
186
187      if ((itage.lt.7200).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
188           (xl.gt.real(numxgrid-1)-0.5).or. &
189           (yl.gt.real(numygrid-1)-0.5)) then             ! no kernel, direct attribution to grid cell
190        if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
191             (jy.le.numygrid-1)) then
192          do ks=1,nspec
193!        print*,'bef ks'
194!         print*,ix,jy,kz,ks,nrelpointer,nclass(i),nage
195          incl=nclass(i)
196             gridunc(ix,jy,kz,ks,nrelpointer,incl,nage)= &
197                  gridunc(ix,jy,kz,ks,nrelpointer,incl,nage)+ &
198                  xmass1(i,ks)/rhoi*weight
199          end do
200        endif
201
202      else                                 ! attribution via uniform kernel
203
204        ddx=xl-real(ix)                   ! distance to left cell border
205        ddy=yl-real(jy)                   ! distance to lower cell border
206        if (ddx.gt.0.5) then
207          ixp=ix+1
208          wx=1.5-ddx
209        else
210          ixp=ix-1
211          wx=0.5+ddx
212        endif
213
214        if (ddy.gt.0.5) then
215          jyp=jy+1
216          wy=1.5-ddy
217        else
218          jyp=jy-1
219          wy=0.5+ddy
220        endif
221
222
223  ! Determine mass fractions for four grid points
224  !**********************************************
225
226        if ((ix.ge.0).and.(ix.le.numxgrid-1)) then
227          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
228            w=wx*wy
229            do ks=1,nspec
230              gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
231                   gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
232                   xmass1(i,ks)/rhoi*weight*w
233            end do
234          endif
235
236          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
237            w=wx*(1.-wy)
238            do ks=1,nspec
239              gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
240                   gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
241                   xmass1(i,ks)/rhoi*weight*w
242            end do
243          endif
244        endif
245
246
247        if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then
248          if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then
249            w=(1.-wx)*(1.-wy)
250            do ks=1,nspec
251              gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
252                   gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
253                   xmass1(i,ks)/rhoi*weight*w
254            end do
255          endif
256
257          if ((jy.ge.0).and.(jy.le.numygrid-1)) then
258            w=(1.-wx)*wy
259            do ks=1,nspec
260              gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
261                   gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
262                   xmass1(i,ks)/rhoi*weight*w
263            end do
264          endif
265        endif
266      endif
267
268
269
270  !************************************
271  ! Do everything for the nested domain
272  !************************************
273
274      if (nested_output.eq.1) then
275        xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
276        yl=(ytra1(i)*dy+youtshiftn)/dyoutn
277        ix=int(xl)
278        if (xl.lt.0.) ix=ix-1
279        jy=int(yl)
280        if (yl.lt.0.) jy=jy-1
281
282
283  ! For particles aged less than 3 hours, attribute particle mass to grid cell
284  ! it resides in rather than use the kernel, in order to avoid its smoothing effect.
285  ! For older particles, use the uniform kernel.
286  ! If a particle is close to the domain boundary, do not use the kernel either.
287  !*****************************************************************************
288
289        if ((itage.lt.7200).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
290             (xl.gt.real(numxgridn-1)-0.5).or. &
291             (yl.gt.real(numygridn-1)-0.5)) then             ! no kernel, direct attribution to grid cell
292          if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
293               (jy.le.numygridn-1)) then
294            do ks=1,nspec
295              griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
296                   griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
297                   xmass1(i,ks)/rhoi*weight
298            end do
299          endif
300
301        else                                 ! attribution via uniform kernel
302
303          ddx=xl-real(ix)                   ! distance to left cell border
304          ddy=yl-real(jy)                   ! distance to lower cell border
305          if (ddx.gt.0.5) then
306            ixp=ix+1
307            wx=1.5-ddx
308          else
309            ixp=ix-1
310            wx=0.5+ddx
311          endif
312
313          if (ddy.gt.0.5) then
314            jyp=jy+1
315            wy=1.5-ddy
316          else
317            jyp=jy-1
318            wy=0.5+ddy
319          endif
320
321
322  ! Determine mass fractions for four grid points
323  !**********************************************
324
325          if ((ix.ge.0).and.(ix.le.numxgridn-1)) then
326            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
327              w=wx*wy
328              do ks=1,nspec
329                griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
330                     griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
331                     xmass1(i,ks)/rhoi*weight*w
332              end do
333            endif
334
335            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
336              w=wx*(1.-wy)
337              do ks=1,nspec
338                griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
339                     griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
340                     xmass1(i,ks)/rhoi*weight*w
341              end do
342            endif
343          endif
344
345
346          if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then
347            if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then
348              w=(1.-wx)*(1.-wy)
349              do ks=1,nspec
350                griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
351                     griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
352                     xmass1(i,ks)/rhoi*weight*w
353              end do
354            endif
355
356            if ((jy.ge.0).and.(jy.le.numygridn-1)) then
357              w=(1.-wx)*wy
358              do ks=1,nspec
359                griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
360                     griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
361                     xmass1(i,ks)/rhoi*weight*w
362              end do
363            endif
364          endif
365        endif
366
367      endif
368    endif
36920  continue
370  end do
371
372  !***********************************************************************
373  ! 2. Evaluate concentrations at receptor points, using the kernel method
374  !***********************************************************************
375
376  do n=1,numreceptor
377
378
379  ! Reset concentrations
380  !*********************
381
382    do ks=1,nspec
383      c(ks)=0.
384    end do
385
386
387  ! Estimate concentration at receptor
388  !***********************************
389
390    do i=1,numpart
391
392      if (itra1(i).ne.itime) goto 40
393      itage=abs(itra1(i)-itramem(i))
394
395      hz=min(50.+0.3*sqrt(real(itage)),hzmax)
396      zd=ztra1(i)/hz
397      if (zd.gt.1.) goto 40          ! save computing time, leave loop
398
399      hx=min((0.29+2.222e-3*sqrt(real(itage)))*dx+ &
400           real(itage)*1.2e-5,hxmax)                     ! 80 km/day
401      xd=(xtra1(i)-xreceptor(n))/hx
402      if (xd*xd.gt.1.) goto 40       ! save computing time, leave loop
403
404      hy=min((0.18+1.389e-3*sqrt(real(itage)))*dy+ &
405           real(itage)*7.5e-6,hymax)                     ! 80 km/day
406      yd=(ytra1(i)-yreceptor(n))/hy
407      if (yd*yd.gt.1.) goto 40       ! save computing time, leave loop
408      h=hx*hy*hz
409
410      r2=xd*xd+yd*yd+zd*zd
411      if (r2.lt.1.) then
412        xkern=factor*(1.-r2)
413        do ks=1,nspec
414          c(ks)=c(ks)+xmass1(i,ks)*xkern/h
415        end do
416      endif
41740    continue
418    end do
419
420    do ks=1,nspec
421      creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n)
422    end do
423  end do
424
425end subroutine conccalc_irreg
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG