source: flexpart.git/src/releaseparticles.f90 @ 7999df47

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 7999df47 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago

Added code, makefile for dev branch

  • Property mode set to 100644
File size: 16.6 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  use random_mod, only: ran1
53
54  implicit none
55
56  !real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint)
57  real :: xaux,yaux,zaux,rfraction
58  real :: topo,rhoaux(2),r,t,rhoout,ddx,ddy,rddx,rddy,p1,p2,p3,p4
59  real :: dz1,dz2,dz,xtn,ytn,xlonav,timecorrect(maxspec),press,pressold
60  real :: presspart,average_timecorrect
61  integer :: itime,numrel,i,j,k,n,ix,jy,ixp,jyp,ipart,minpart,ii
62  integer :: indz,indzp,kz,ngrid
63  integer :: nweeks,ndayofweek,nhour,jjjjmmdd,ihmmss,mm
64  real(kind=dp) :: juldate,julmonday,jul,jullocal,juldiff
65  real,parameter :: eps=nxmax/3.e5,eps2=1.e-6
66
67  integer :: idummy = -7
68  !save idummy,xmasssave
69  !data idummy/-7/,xmasssave/maxpoint*0./
70
71
72
73  ! Determine the actual date and time in Greenwich (i.e., UTC + correction for daylight savings time)
74  !*****************************************************************************
75
76  julmonday=juldate(19000101,0)          ! this is a Monday
77  jul=bdate+real(itime,kind=dp)/86400._dp    ! this is the current day
78  call caldate(jul,jjjjmmdd,ihmmss)
79  mm=(jjjjmmdd-10000*(jjjjmmdd/10000))/100
80  if ((mm.ge.4).and.(mm.le.9)) jul=jul+1._dp/24._dp   ! daylight savings time in summer
81
82
83  ! For every release point, check whether we are in the release time interval
84  !***************************************************************************
85
86  minpart=1
87  do i=1,numpoint
88    if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
89         (itime.le.ireleaseend(i))) then
90
91  ! Determine the local day and time
92  !*********************************
93
94      xlonav=xlon0+(xpoint2(i)+xpoint1(i))/2.*dx  ! longitude needed to determine local time
95      if (xlonav.lt.-180.) xlonav=xlonav+360.
96      if (xlonav.gt.180.) xlonav=xlonav-360.
97      jullocal=jul+real(xlonav,kind=dp)/360._dp   ! correct approximately for time zone to obtain local time
98
99      juldiff=jullocal-julmonday
100      nweeks=int(juldiff/7._dp)
101      juldiff=juldiff-real(nweeks,kind=dp)*7._dp
102      ndayofweek=int(juldiff)+1              ! this is the current day of week, starting with Monday
103      nhour=nint((juldiff-real(ndayofweek-1,kind=dp))*24._dp)    ! this is the current hour
104      if (nhour.eq.0) then
105        nhour=24
106        ndayofweek=ndayofweek-1
107        if (ndayofweek.eq.0) ndayofweek=7
108      endif
109
110  ! Calculate a species- and time-dependent correction factor, distinguishing between
111  ! area (those with release starting at surface) and point (release starting above surface) sources
112  ! Also, calculate an average time correction factor (species independent)
113  !*****************************************************************************
114      average_timecorrect=0.
115      do k=1,nspec
116        if (zpoint1(i).gt.0.5) then      ! point source
117          timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
118        else                             ! area source
119          timecorrect(k)=area_hour(k,nhour)*area_dow(k,ndayofweek)
120        endif
121        average_timecorrect=average_timecorrect+timecorrect(k)
122      end do
123      average_timecorrect=average_timecorrect/real(nspec)
124
125  ! Determine number of particles to be released this time; at start and at end of release,
126  ! only half the particles are released
127  !*****************************************************************************
128
129      if (ireleasestart(i).ne.ireleaseend(i)) then
130        rfraction=abs(real(npart(i))*real(lsynctime)/ &
131             real(ireleaseend(i)-ireleasestart(i)))
132        if ((itime.eq.ireleasestart(i)).or. &
133             (itime.eq.ireleaseend(i))) rfraction=rfraction/2.
134
135  ! Take the species-average time correction factor in order to scale the
136  ! number of particles released this time
137  !**********************************************************************
138        rfraction=rfraction*average_timecorrect
139
140        rfraction=rfraction+xmasssave(i)  ! number to be released at this time
141        numrel=int(rfraction)
142        xmasssave(i)=rfraction-real(numrel)
143      else
144        numrel=npart(i)
145      endif
146
147      xaux=xpoint2(i)-xpoint1(i)
148      yaux=ypoint2(i)-ypoint1(i)
149      zaux=zpoint2(i)-zpoint1(i)
150      do j=1,numrel                       ! loop over particles to be released this time
151        do ipart=minpart,maxpart          ! search for free storage space
152
153  ! If a free storage space is found, attribute everything to this array element
154  !*****************************************************************************
155
156          if (itra1(ipart).ne.itime) then
157
158  ! Particle coordinates are determined by using a random position within the release volume
159  !*****************************************************************************
160
161  ! Determine horizontal particle position
162  !***************************************
163
164            xtra1(ipart)=xpoint1(i)+ran1(idummy)*xaux
165            if (xglobal) then
166              if (xtra1(ipart).gt.real(nxmin1)) xtra1(ipart)= &
167                   xtra1(ipart)-real(nxmin1)
168              if (xtra1(ipart).lt.0.) xtra1(ipart)= &
169                   xtra1(ipart)+real(nxmin1)
170            endif
171            ytra1(ipart)=ypoint1(i)+ran1(idummy)*yaux
172
173  ! Assign mass to particle: Total mass divided by total number of particles.
174  ! Time variation has partly been taken into account already by a species-average
175  ! correction factor, by which the number of particles released this time has been
176  ! scaled. Adjust the mass per particle by the species-dependent time correction factor
177  ! divided by the species-average one
178  !*****************************************************************************
179            do k=1,nspec
180               xmass1(ipart,k)=xmass(i,k)/real(npart(i)) &
181                    *timecorrect(k)/average_timecorrect
182  !            write (*,*) 'xmass1: ',xmass1(ipart,k),ipart,k
183  ! Assign certain properties to particle
184  !**************************************
185            end do
186            nclass(ipart)=min(int(ran1(idummy)*real(nclassunc))+1, &
187                 nclassunc)
188            numparticlecount=numparticlecount+1
189            if (mquasilag.eq.0) then
190              npoint(ipart)=i
191            else
192              npoint(ipart)=numparticlecount
193            endif
194            idt(ipart)=mintime               ! first time step
195            itra1(ipart)=itime
196            itramem(ipart)=itra1(ipart)
197            itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
198
199
200  ! Determine vertical particle position
201  !*************************************
202
203            ztra1(ipart)=zpoint1(i)+ran1(idummy)*zaux
204
205  ! Interpolation of topography and density
206  !****************************************
207
208  ! Determine the nest we are in
209  !*****************************
210
211            ngrid=0
212            do k=numbnests,1,-1
213              if ((xtra1(ipart).gt.xln(k)+eps).and. &
214                   (xtra1(ipart).lt.xrn(k)-eps).and. &
215                   (ytra1(ipart).gt.yln(k)+eps).and. &
216                   (ytra1(ipart).lt.yrn(k)-eps)) then
217                ngrid=k
218                goto 43
219              endif
220            end do
22143          continue
222
223  ! Determine (nested) grid coordinates and auxiliary parameters used for interpolation
224  !*****************************************************************************
225
226            if (ngrid.gt.0) then
227              xtn=(xtra1(ipart)-xln(ngrid))*xresoln(ngrid)
228              ytn=(ytra1(ipart)-yln(ngrid))*yresoln(ngrid)
229              ix=int(xtn)
230              jy=int(ytn)
231              ddy=ytn-real(jy)
232              ddx=xtn-real(ix)
233            else
234              ix=int(xtra1(ipart))
235              jy=int(ytra1(ipart))
236              ddy=ytra1(ipart)-real(jy)
237              ddx=xtra1(ipart)-real(ix)
238            endif
239            ixp=ix+1
240            jyp=jy+1
241            rddx=1.-ddx
242            rddy=1.-ddy
243            p1=rddx*rddy
244            p2=ddx*rddy
245            p3=rddx*ddy
246            p4=ddx*ddy
247
248            if (ngrid.gt.0) then
249              topo=p1*oron(ix ,jy ,ngrid) &
250                   + p2*oron(ixp,jy ,ngrid) &
251                   + p3*oron(ix ,jyp,ngrid) &
252                   + p4*oron(ixp,jyp,ngrid)
253            else
254              topo=p1*oro(ix ,jy) &
255                   + p2*oro(ixp,jy) &
256                   + p3*oro(ix ,jyp) &
257                   + p4*oro(ixp,jyp)
258            endif
259
260  ! If starting height is in pressure coordinates, retrieve pressure profile and convert zpart1 to meters
261  !*****************************************************************************
262            if (kindz(i).eq.3) then
263              presspart=ztra1(ipart)
264              do kz=1,nz
265                if (ngrid.gt.0) then
266                  r=p1*rhon(ix ,jy ,kz,2,ngrid) &
267                       +p2*rhon(ixp,jy ,kz,2,ngrid) &
268                       +p3*rhon(ix ,jyp,kz,2,ngrid) &
269                       +p4*rhon(ixp,jyp,kz,2,ngrid)
270                  t=p1*ttn(ix ,jy ,kz,2,ngrid) &
271                       +p2*ttn(ixp,jy ,kz,2,ngrid) &
272                       +p3*ttn(ix ,jyp,kz,2,ngrid) &
273                       +p4*ttn(ixp,jyp,kz,2,ngrid)
274                else
275                  r=p1*rho(ix ,jy ,kz,2) &
276                       +p2*rho(ixp,jy ,kz,2) &
277                       +p3*rho(ix ,jyp,kz,2) &
278                       +p4*rho(ixp,jyp,kz,2)
279                  t=p1*tt(ix ,jy ,kz,2) &
280                       +p2*tt(ixp,jy ,kz,2) &
281                       +p3*tt(ix ,jyp,kz,2) &
282                       +p4*tt(ixp,jyp,kz,2)
283                endif
284                press=r*r_air*t/100.
285                if (kz.eq.1) pressold=press
286
287                if (press.lt.presspart) then
288                  if (kz.eq.1) then
289                    ztra1(ipart)=height(1)/2.
290                  else
291                    dz1=pressold-presspart
292                    dz2=presspart-press
293                    ztra1(ipart)=(height(kz-1)*dz2+height(kz)*dz1) &
294                         /(dz1+dz2)
295                  endif
296                  goto 71
297                endif
298                pressold=press
299              end do
30071            continue
301            endif
302
303  ! If release positions are given in meters above sea level, subtract the
304  ! topography from the starting height
305  !***********************************************************************
306
307            if (kindz(i).eq.2) ztra1(ipart)=ztra1(ipart)-topo
308            if (ztra1(ipart).lt.eps2) ztra1(ipart)=eps2   ! Minimum starting height is eps2
309            if (ztra1(ipart).gt.height(nz)-0.5) ztra1(ipart)= &
310                 height(nz)-0.5 ! Maximum starting height is uppermost level - 0.5 meters
311
312
313
314  ! For special simulations, multiply particle concentration air density;
315  ! Simply take the 2nd field in memory to do this (accurate enough)
316  !***********************************************************************
317  !AF IND_SOURCE switches between different units for concentrations at the source
318  !Af    NOTE that in backward simulations the release of particles takes place at the
319  !Af         receptor and the sampling at the source.
320  !Af          1="mass"
321  !Af          2="mass mixing ratio"
322  !Af IND_RECEPTOR switches between different units for concentrations at the receptor
323  !Af          1="mass"
324  !Af          2="mass mixing ratio"
325
326  !Af switches for the releasefile:
327  !Af IND_REL =  1 : xmass * rho
328  !Af IND_REL =  0 : xmass * 1
329
330  !Af ind_rel is defined in readcommand.f
331
332            if (ind_rel .eq. 1) then
333
334  ! Interpolate the air density
335  !****************************
336
337              do ii=2,nz
338                if (height(ii).gt.ztra1(ipart)) then
339                  indz=ii-1
340                  indzp=ii
341                  goto 6
342                endif
343              end do
3446             continue
345
346              dz1=ztra1(ipart)-height(indz)
347              dz2=height(indzp)-ztra1(ipart)
348              dz=1./(dz1+dz2)
349
350              if (ngrid.gt.0) then
351                do n=1,2
352                  rhoaux(n)=p1*rhon(ix ,jy ,indz+n-1,2,ngrid) &
353                       +p2*rhon(ixp,jy ,indz+n-1,2,ngrid) &
354                       +p3*rhon(ix ,jyp,indz+n-1,2,ngrid) &
355                       +p4*rhon(ixp,jyp,indz+n-1,2,ngrid)
356                end do
357              else
358                do n=1,2
359                  rhoaux(n)=p1*rho(ix ,jy ,indz+n-1,2) &
360                       +p2*rho(ixp,jy ,indz+n-1,2) &
361                       +p3*rho(ix ,jyp,indz+n-1,2) &
362                       +p4*rho(ixp,jyp,indz+n-1,2)
363                end do
364              endif
365              rhoout=(dz2*rhoaux(1)+dz1*rhoaux(2))*dz
366              rho_rel(i)=rhoout
367
368
369  ! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density
370  !********************************************************************
371
372              do k=1,nspec
373                xmass1(ipart,k)=xmass1(ipart,k)*rhoout
374              end do
375            endif
376
377
378            numpart=max(numpart,ipart)
379            goto 34      ! Storage space has been found, stop searching
380          endif
381        end do
382        if (ipart.gt.maxpart) goto 996
383
38434      minpart=ipart+1
385      end do
386      endif
387  end do
388
389
390  return
391
392996   continue
393  write(*,*) '#####################################################'
394  write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####'
395  write(*,*) '####                                             ####'
396  write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED  ####'
397  write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE  ####'
398  write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
399  write(*,*) '#### OR REDUCE NUMBER OF RELEASE POINTS.         ####'
400  write(*,*) '#####################################################'
401  stop
402
403end subroutine releaseparticles
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG