source: trunk/src/releaseparticles.f90 @ 4

Last change on this file since 4 was 4, checked in by mlanger, 10 years ago
File size: 16.5 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine releaseparticles(itime)
23  !                              o
24  !*****************************************************************************
25  !                                                                            *
26  !     This subroutine releases particles from the release locations.         *
27  !                                                                            *
28  !     It searches for a "vacant" storage space and assigns all particle      *
29  !     information to that space. A space is vacant either when no particle   *
30  !     is yet assigned to it, or when it's particle is expired and, thus,     *
31  !     the storage space is made available to a new particle.                 *
32  !                                                                            *
33  !     Author: A. Stohl                                                       *
34  !                                                                            *
35  !     29 June 2002                                                           *
36  !                                                                            *
37  !*****************************************************************************
38  !                                                                            *
39  ! Variables:                                                                 *
40  ! itime [s]            current time                                          *
41  ! ireleasestart, ireleaseend          start and end times of all releases    *
42  ! npart(maxpoint)      number of particles to be released in total           *
43  ! numrel               number of particles to be released during this time   *
44  !                      step                                                  *
45  !                                                                            *
46  !*****************************************************************************
47
48  use point_mod
49  use xmass_mod
50  use par_mod
51  use com_mod
52
53  implicit none
54
55  !real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint)
56  real :: xaux,yaux,zaux,ran1,rfraction
57  real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4
58  real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold
59  real :: presspart,average_timecorrect
60  integer :: itime,numrel,i,j,k,n,ix,jy,ixp,jyp,ipart,minpart,ii
61  integer :: indz,indzp,kz,ngrid
62  integer :: nweeks,ndayofweek,nhour,jjjjmmdd,ihmmss,mm
63  real(kind=dp) :: juldate,julmonday,jul,jullocal,juldiff
64  real,parameter :: eps=nxmax/3.e5,eps2=1.e-6
65
66  integer :: idummy = -7
67  !save idummy,xmasssave
68  !data idummy/-7/,xmasssave/maxpoint*0./
69
70
71  ! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time)
72  !*****************************************************************************
73
74  julmonday=juldate(19000101,0)          ! this is a Monday
75  jul=bdate+real(itime,kind=dp)/86400._dp    ! this is the current day
76  call caldate(jul,jjjjmmdd,ihmmss)
77  mm=(jjjjmmdd-10000*(jjjjmmdd/10000))/100
78  if ((mm.ge.4).and.(mm.le.9)) jul=jul+1._dp/24._dp   ! daylight savings time in summer
79
80
81  ! For every release point, check whether we are in the release time interval
82  !***************************************************************************
83
84  minpart=1
85  do i=1,numpoint
86    if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
87         (itime.le.ireleaseend(i))) then
88
89  ! Determine the local day and time
90  !*********************************
91
92      xlonav=xlon0+(xpoint2(i)+xpoint1(i))/2.*dx  ! longitude needed to determine local time
93      if (xlonav.lt.-180.) xlonav=xlonav+360.
94      if (xlonav.gt.180.) xlonav=xlonav-360.
95      jullocal=jul+real(xlonav,kind=dp)/360._dp   ! correct approximately for time zone to obtain local time
96
97      juldiff=jullocal-julmonday
98      nweeks=int(juldiff/7._dp)
99      juldiff=juldiff-real(nweeks,kind=dp)*7._dp
100      ndayofweek=int(juldiff)+1              ! this is the current day of week, starting with Monday
101      nhour=nint((juldiff-real(ndayofweek-1,kind=dp))*24._dp)    ! this is the current hour
102      if (nhour.eq.0) then
103        nhour=24
104        ndayofweek=ndayofweek-1
105        if (ndayofweek.eq.0) ndayofweek=7
106      endif
107
108  ! Calculate a species- and time-dependent correction factor, distinguishing between
109  ! area (those with release starting at surface) and point (release starting above surface) sources
110  ! Also, calculate an average time correction factor (species independent)
111  !*****************************************************************************
112      average_timecorrect=0.
113      do k=1,nspec
114        if (zpoint1(i).gt.0.5) then      ! point source
115          timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
116        else                             ! area source
117          timecorrect(k)=area_hour(k,nhour)*area_dow(k,ndayofweek)
118        endif
119        average_timecorrect=average_timecorrect+timecorrect(k)
120      end do
121      average_timecorrect=average_timecorrect/real(nspec)
122
123  ! Determine number of particles to be released this time; at start and at end of release,
124  ! only half the particles are released
125  !*****************************************************************************
126
127      if (ireleasestart(i).ne.ireleaseend(i)) then
128        rfraction=abs(real(npart(i))*real(lsynctime)/ &
129             real(ireleaseend(i)-ireleasestart(i)))
130        if ((itime.eq.ireleasestart(i)).or. &
131             (itime.eq.ireleaseend(i))) rfraction=rfraction/2.
132
133  ! Take the species-average time correction factor in order to scale the
134  ! number of particles released this time
135  !**********************************************************************
136        rfraction=rfraction*average_timecorrect
137
138        rfraction=rfraction+xmasssave(i)  ! number to be released at this time
139        numrel=int(rfraction)
140        xmasssave(i)=rfraction-real(numrel)
141      else
142        numrel=npart(i)
143      endif
144
145      xaux=xpoint2(i)-xpoint1(i)
146      yaux=ypoint2(i)-ypoint1(i)
147      zaux=zpoint2(i)-zpoint1(i)
148      do j=1,numrel                       ! loop over particles to be released this time
149        do ipart=minpart,maxpart          ! search for free storage space
150
151  ! If a free storage space is found, attribute everything to this array element
152  !*****************************************************************************
153
154          if (itra1(ipart).ne.itime) then
155
156  ! Particle coordinates are determined by using a random position within the release volume
157  !*****************************************************************************
158
159  ! Determine horizontal particle position
160  !***************************************
161
162            xtra1(ipart)=xpoint1(i)+ran1(idummy)*xaux
163            if (xglobal) then
164              if (xtra1(ipart).gt.real(nxmin1)) xtra1(ipart)= &
165                   xtra1(ipart)-real(nxmin1)
166              if (xtra1(ipart).lt.0.) xtra1(ipart)= &
167                   xtra1(ipart)+real(nxmin1)
168            endif
169            ytra1(ipart)=ypoint1(i)+ran1(idummy)*yaux
170
171  ! Assign mass to particle: Total mass divided by total number of particles.
172  ! Time variation has partly been taken into account already by a species-average
173  ! correction factor, by which the number of particles released this time has been
174  ! scaled. Adjust the mass per particle by the species-dependent time correction factor
175  ! divided by the species-average one
176  !*****************************************************************************
177            do k=1,nspec
178               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
179                    *timecorrect(k)/average_timecorrect
180  !            write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
181  ! Assign certain properties to particle
182  !**************************************
183            end do
184            nclass(ipart)=min(int(ran1(idummy)*real(nclassunc))+1, &
185                 nclassunc)
186            numparticlecount=numparticlecount+1
187            if (mquasilag.eq.0) then
188              npoint(ipart)=i
189            else
190              npoint(ipart)=numparticlecount
191            endif
192            idt(ipart)=mintime               ! first time step
193            itra1(ipart)=itime
194            itramem(ipart)=itra1(ipart)
195            itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
196
197
198  ! Determine vertical particle position
199  !*************************************
200
201            ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
202
203  ! Interpolation of topography and density
204  !****************************************
205
206  ! Determine the nest we are in
207  !*****************************
208
209            ngrid=0
210            do k=numbnests,1,-1
211              if ((xtra1(ipart).gt.xln(k)+eps).and. &
212                   (xtra1(ipart).lt.xrn(k)-eps).and. &
213                   (ytra1(ipart).gt.yln(k)+eps).and. &
214                   (ytra1(ipart).lt.yrn(k)-eps)) then
215                ngrid=k
216                goto 43
217              endif
218            end do
21943          continue
220
221  ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
222  !*****************************************************************************
223
224            if (ngrid.gt.0) then
225              xtn=(xtra1(ipart)-xln(ngrid))*xresoln(ngrid)
226              ytn=(ytra1(ipart)-yln(ngrid))*yresoln(ngrid)
227              ix=int(xtn)
228              jy=int(ytn)
229              ddy=ytn-real(jy)
230              ddx=xtn-real(ix)
231            else
232              ix=int(xtra1(ipart))
233              jy=int(ytra1(ipart))
234              ddy=ytra1(ipart)-real(jy)
235              ddx=xtra1(ipart)-real(ix)
236            endif
237            ixp=ix+1
238            jyp=jy+1
239            rddx=1.-ddx
240            rddy=1.-ddy
241            p1=rddx*rddy
242            p2=ddx*rddy
243            p3=rddx*ddy
244            p4=ddx*ddy
245
246            if (ngrid.gt.0) then
247              topo=p1*oron(ix ,jy ,ngrid) &
248                   + p2*oron(ixp,jy ,ngrid) &
249                   + p3*oron(ix ,jyp,ngrid) &
250                   + p4*oron(ixp,jyp,ngrid)
251            else
252              topo=p1*oro(ix ,jy) &
253                   + p2*oro(ixp,jy) &
254                   + p3*oro(ix ,jyp) &
255                   + p4*oro(ixp,jyp)
256            endif
257
258  ! If starting height is in pressure coordinates, retrieve pressure profile and convert zpart1 to meters
259  !*****************************************************************************
260            if (kindz(i).eq.3) then
261              presspart=ztra1(ipart)
262              do kz=1,nz
263                if (ngrid.gt.0) then
264                  r=p1*rhon(ix ,jy ,kz,2,ngrid) &
265                       +p2*rhon(ixp,jy ,kz,2,ngrid) &
266                       +p3*rhon(ix ,jyp,kz,2,ngrid) &
267                       +p4*rhon(ixp,jyp,kz,2,ngrid)
268                  t=p1*ttn(ix ,jy ,kz,2,ngrid) &
269                       +p2*ttn(ixp,jy ,kz,2,ngrid) &
270                       +p3*ttn(ix ,jyp,kz,2,ngrid) &
271                       +p4*ttn(ixp,jyp,kz,2,ngrid)
272                else
273                  r=p1*rho(ix ,jy ,kz,2) &
274                       +p2*rho(ixp,jy ,kz,2) &
275                       +p3*rho(ix ,jyp,kz,2) &
276                       +p4*rho(ixp,jyp,kz,2)
277                  t=p1*tt(ix ,jy ,kz,2) &
278                       +p2*tt(ixp,jy ,kz,2) &
279                       +p3*tt(ix ,jyp,kz,2) &
280                       +p4*tt(ixp,jyp,kz,2)
281                endif
282                press=r*r_air*t/100.
283                if (kz.eq.1) pressold=press
284
285                if (press.lt.presspart) then
286                  if (kz.eq.1) then
287                    ztra1(ipart)=height(1)/2.
288                  else
289                    dz1=pressold-presspart
290                    dz2=presspart-press
291                    ztra1(ipart)=(height(kz-1)*dz2+height(kz)*dz1) &
292                         /(dz1+dz2)
293                  endif
294                  goto 71
295                endif
296                pressold=press
297              end do
29871            continue
299            endif
300
301  ! If release positions are given in meters above sea level, subtract the
302  ! topography from the starting height
303  !***********************************************************************
304
305            if (kindz(i).eq.2) ztra1(ipart)=ztra1(ipart)-topo
306            if (ztra1(ipart).lt.eps2) ztra1(ipart)=eps2   ! Minimum starting height is eps2
307            if (ztra1(ipart).gt.height(nz)-0.5) ztra1(ipart)= &
308                 height(nz)-0.5 ! Maximum starting height is uppermost level - 0.5 meters
309
310
311
312  ! For special simulations, multiply particle concentration air density;
313  ! Simply take the 2nd field in memory to do this (accurate enough)
314  !***********************************************************************
315  !AF IND_SOURCE switches between different units for concentrations at the source
316  !Af    NOTE that in backward simulations the release of particles takes place at the
317  !Af         receptor and the sampling at the source.
318  !Af          1="mass"
319  !Af          2="mass mixing ratio"
320  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
321  !Af          1="mass"
322  !Af          2="mass mixing ratio"
323
324  !Af switches for the releasefile:
325  !Af IND_REL =  1 : xmass * rho
326  !Af IND_REL =  0 : xmass * 1
327
328  !Af ind_rel is defined in readcommand.f
329
330            if (ind_rel .eq. 1) then
331
332  ! Interpolate the air density
333  !****************************
334
335              do ii=2,nz
336                if (height(ii).gt.ztra1(ipart)) then
337                  indz=ii-1
338                  indzp=ii
339                  goto 6
340                endif
341              end do
3426             continue
343
344              dz1=ztra1(ipart)-height(indz)
345              dz2=height(indzp)-ztra1(ipart)
346              dz=1./(dz1+dz2)
347
348              if (ngrid.gt.0) then
349                do n=1,2
350                  rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,2,ngrid) &
351                       +p2*rhon(ixp,jy ,indz+n-1,2,ngrid) &
352                       +p3*rhon(ix ,jyp,indz+n-1,2,ngrid) &
353                       +p4*rhon(ixp,jyp,indz+n-1,2,ngrid)
354                end do
355              else
356                do n=1,2
357                  rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,2) &
358                       +p2*rho(ixp,jy ,indz+n-1,2) &
359                       +p3*rho(ix ,jyp,indz+n-1,2) &
360                       +p4*rho(ixp,jyp,indz+n-1,2)
361                end do
362              endif
363              rhoout=(dz2*rhoaux(1)+dz1*rhoaux(2))*dz
364              rho_rel(i)=rhoout
365
366
367  ! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density
368  !********************************************************************
369
370              do k=1,nspec
371                xmass1(ipart,k)=xmass1(ipart,k)*rhoout
372              end do
373            endif
374
375
376            numpart=max(numpart,ipart)
377            goto 34      ! Storage space has been found, stop searching
378          endif
379        end do
380        if (ipart.gt.maxpart) goto 996
381
38234      minpart=ipart+1
383      end do
384      endif
385  end do
386
387
388  return
389
390996   continue
391  write(*,*) '#####################################################'
392  write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####'
393  write(*,*) '####                                             ####'
394  write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED  ####'
395  write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE  ####'
396  write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
397  write(*,*) '#### OR REDUCE NUMBER OF RELEASE POINTS.         ####'
398  write(*,*) '#####################################################'
399  stop
400
401end subroutine releaseparticles
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG