source: flexpart.git/src/releaseparticles_mpi.f90 @ 16b61a5

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 16b61a5 was 16b61a5, checked in by Espen Sollum ATMOS <eso@…>, 7 years ago

Reworked the domain-filling option (MPI). Fixed a slow loop which had errors in loop counter (MPI)

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