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
66 integer :: idummy = -7
75 jul=bdate+
real(itime,kind=dp)/86400._dp
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
86 if ((itime.ge.ireleasestart(i)).and. &
87 (itime.le.ireleaseend(i)))
then
92 xlonav=xlon0+(xpoint2(i)+xpoint1(i))/2.*dx
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
97 juldiff=jullocal-julmonday
98 nweeks=int(juldiff/7._dp)
99 juldiff=juldiff-
real(nweeks,kind=dp)*7._dp
100 ndayofweek=int(juldiff)+1
101 nhour=nint((juldiff-
real(ndayofweek-1,kind=dp))*24._dp)
104 ndayofweek=ndayofweek-1
105 if (ndayofweek.eq.0) ndayofweek=7
112 average_timecorrect=0.
114 if (zpoint1(i).gt.0.5)
then
115 timecorrect(k)=point_hour(k,nhour)*point_dow(k,ndayofweek)
117 timecorrect(k)=area_hour(k,nhour)*area_dow(k,ndayofweek)
119 average_timecorrect=average_timecorrect+timecorrect(k)
121 average_timecorrect=average_timecorrect/
real(nspec)
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.
136 rfraction=rfraction*average_timecorrect
138 rfraction=rfraction+xmasssave(i)
139 numrel=int(rfraction)
140 xmasssave(i)=rfraction-
real(numrel)
145 xaux=xpoint2(i)-xpoint1(i)
146 yaux=ypoint2(i)-ypoint1(i)
147 zaux=zpoint2(i)-zpoint1(i)
149 do ipart=minpart,maxpart
154 if (itra1(ipart).ne.itime)
then
162 xtra1(ipart)=xpoint1(i)+
ran1(idummy)*xaux
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)
169 ytra1(ipart)=ypoint1(i)+
ran1(idummy)*yaux
178 xmass1(ipart,k)=xmass(i,k)/
real(npart(i)) &
179 *timecorrect(k)/average_timecorrect
184 nclass(ipart)=min(int(
ran1(idummy)*
real(nclassunc))+1, &
186 numparticlecount=numparticlecount+1
187 if (mquasilag.eq.0)
then
190 npoint(ipart)=numparticlecount
194 itramem(ipart)=itra1(ipart)
195 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
201 ztra1(ipart)=zpoint1(i)+
ran1(idummy)*zaux
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
225 xtn=(xtra1(ipart)-xln(ngrid))*xresoln(ngrid)
226 ytn=(ytra1(ipart)-yln(ngrid))*yresoln(ngrid)
234 ddy=ytra1(ipart)-
real(jy)
235 ddx=xtra1(ipart)-
real(ix)
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)
252 topo=p1*oro(ix ,jy) &
260 if (kindz(i).eq.3)
then
261 presspart=ztra1(ipart)
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)
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) &
283 if (kz.eq.1) pressold=press
285 if (press.lt.presspart)
then
287 ztra1(ipart)=height(1)/2.
289 dz1=pressold-presspart
291 ztra1(ipart)=(height(kz-1)*dz2+height(kz)*dz1) &
305 if (kindz(i).eq.2) ztra1(ipart)=ztra1(ipart)-topo
306 if (ztra1(ipart).lt.eps2) ztra1(ipart)=eps2
307 if (ztra1(ipart).gt.height(nz)-0.5) ztra1(ipart)= &
330 if (ind_rel .eq. 1)
then
336 if (height(ii).gt.ztra1(ipart))
then
344 dz1=ztra1(ipart)-height(indz)
345 dz2=height(indzp)-ztra1(ipart)
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)
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)
363 rhoout=(dz2*rhoaux(1)+dz1*rhoaux(2))*dz
371 xmass1(ipart,k)=xmass1(ipart,k)*rhoout
376 numpart=max(numpart,ipart)
380 if (ipart.gt.maxpart) goto 996
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(*,*)
'#####################################################'
subroutine releaseparticles(itime)
subroutine caldate(juldate, yyyymmdd, hhmiss)
real(kind=dp) function juldate(yyyymmdd, hhmiss)