source: flexpart.git/src/releaseparticles.f90 @ 5f2e8f6

flexpart-noresm
Last change on this file since 5f2e8f6 was 5f2e8f6, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 8 years ago

new flexpart noresm code in src

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