source: flexpart.git/src/releaseparticles.f90 @ 462f74b

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 462f74b was 462f74b, checked in by Sabine Eckhardt <sabine@…>, 8 years ago

first version of the backward scavenging

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