54 integer :: itime,itage,i,ix,jy,ixp,jyp,kz,ks,n,nage
55 integer :: il,ind,indz,indzp,nrelpointer
56 real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz
57 real :: weight,hx,hy,hz,h,xd,yd,zd,xkern,r2,c(maxspec),ddx,ddy
58 real :: rhoprof(2),rhoi
60 real,
parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150.
70 if (itra1(i).ne.itime) goto 20
73 itage=abs(itra1(i)-itramem(i))
75 if (itage.lt.lage(nage)) goto 33
98 if ( ind_samp .eq. -1 )
then
104 ddx=xtra1(i)-
real(ix)
105 ddy=ytra1(i)-
real(jy)
114 if (height(il).gt.ztra1(i))
then
122 dz1=ztra1(i)-height(indz)
123 dz2=height(indzp)-ztra1(i)
129 rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) &
130 +p2*rho(ixp,jy ,ind,2) &
131 +p3*rho(ix ,jyp,ind,2) &
132 +p4*rho(ixp,jyp,ind,2)
134 rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz
135 elseif (ind_samp.eq.0)
then
151 if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1))
then
154 nrelpointer=npoint(i)
158 if (outheight(kz).gt.ztra1(i)) goto 21
161 if (kz.le.numzgrid)
then
168 xl=(xtra1(i)*dx+xoutshift)/dxout
169 yl=(ytra1(i)*dy+youtshift)/dyout
171 if (xl.lt.0.) ix=ix-1
173 if (yl.lt.0.) jy=jy-1
184 if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
185 (xl.gt.
real(numxgrid-1)-0.5).or. &
186 (yl.gt.
real(numygrid-1)-0.5)) then
187 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. &
188 (jy.le.numygrid-1))
then
190 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
191 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
192 xmass1(i,ks)/rhoi*weight
220 if ((ix.ge.0).and.(ix.le.numxgrid-1))
then
221 if ((jy.ge.0).and.(jy.le.numygrid-1))
then
224 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
225 gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
226 xmass1(i,ks)/rhoi*weight*w
230 if ((jyp.ge.0).and.(jyp.le.numygrid-1))
then
233 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
234 gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
235 xmass1(i,ks)/rhoi*weight*w
241 if ((ixp.ge.0).and.(ixp.le.numxgrid-1))
then
242 if ((jyp.ge.0).and.(jyp.le.numygrid-1))
then
245 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
246 gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
247 xmass1(i,ks)/rhoi*weight*w
251 if ((jy.ge.0).and.(jy.le.numygrid-1))
then
254 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
255 gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
256 xmass1(i,ks)/rhoi*weight*w
268 if (nested_output.eq.1)
then
269 xl=(xtra1(i)*dx+xoutshiftn)/dxoutn
270 yl=(ytra1(i)*dy+youtshiftn)/dyoutn
272 if (xl.lt.0.) ix=ix-1
274 if (yl.lt.0.) jy=jy-1
283 if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. &
284 (xl.gt.
real(numxgridn-1)-0.5).or. &
285 (yl.gt.
real(numygridn-1)-0.5)) then
286 if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. &
287 (jy.le.numygridn-1))
then
289 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
290 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
291 xmass1(i,ks)/rhoi*weight
319 if ((ix.ge.0).and.(ix.le.numxgridn-1))
then
320 if ((jy.ge.0).and.(jy.le.numygridn-1))
then
323 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= &
324 griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
325 xmass1(i,ks)/rhoi*weight*w
329 if ((jyp.ge.0).and.(jyp.le.numygridn-1))
then
332 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
333 griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
334 xmass1(i,ks)/rhoi*weight*w
340 if ((ixp.ge.0).and.(ixp.le.numxgridn-1))
then
341 if ((jyp.ge.0).and.(jyp.le.numygridn-1))
then
344 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= &
345 griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ &
346 xmass1(i,ks)/rhoi*weight*w
350 if ((jy.ge.0).and.(jy.le.numygridn-1))
then
353 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= &
354 griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ &
355 xmass1(i,ks)/rhoi*weight*w
386 if (itra1(i).ne.itime) goto 40
387 itage=abs(itra1(i)-itramem(i))
389 hz=min(50.+0.3*sqrt(
real(itage)),hzmax)
391 if (zd.gt.1.) goto 40
393 hx=min((0.29+2.222e-3*sqrt(
real(itage)))*dx+ &
394 real(itage)*1.2e-5,hxmax)
395 xd=(xtra1(i)-xreceptor(n))/hx
396 if (xd*xd.gt.1.) goto 40
398 hy=min((0.18+1.389e-3*sqrt(
real(itage)))*dy+ &
399 real(itage)*7.5e-6,hymax)
400 yd=(ytra1(i)-yreceptor(n))/hy
401 if (yd*yd.gt.1.) goto 40
408 c(ks)=c(ks)+xmass1(i,ks)*xkern/h
415 creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n)
subroutine conccalc(itime, weight)