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

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

sources for flexwrf v3.1

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