!********************************************************************** ! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 * ! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, * ! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann * ! * ! This file is part of FLEXPART. * ! * ! FLEXPART is free software: you can redistribute it and/or modify * ! it under the terms of the GNU General Public License as published by* ! the Free Software Foundation, either version 3 of the License, or * ! (at your option) any later version. * ! * ! FLEXPART is distributed in the hope that it will be useful, * ! but WITHOUT ANY WARRANTY; without even the implied warranty of * ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * ! GNU General Public License for more details. * ! * ! You should have received a copy of the GNU General Public License * ! along with FLEXPART. If not, see . * !********************************************************************** subroutine conccalc(itime,weight) ! i i !***************************************************************************** ! * ! Calculation of the concentrations on a regular grid using volume * ! sampling * ! * ! Author: A. Stohl * ! * ! 24 May 1996 * ! * ! April 2000: Update to calculate age spectra * ! Bug fix to avoid negative conc. at the domain boundaries, * ! as suggested by Petra Seibert * ! * ! 2 July 2002: re-order if-statements in order to optimize CPU time * ! * ! * !***************************************************************************** ! * ! Variables: * ! nspeciesdim = nspec for forward runs, 1 for backward runs * ! * !***************************************************************************** use unc_mod use outg_mod use par_mod use com_mod implicit none integer :: itime,itage,i,ix,jy,ixp,jyp,kz,ks,n,nage integer :: il,ind,indz,indzp,nrelpointer real :: rddx,rddy,p1,p2,p3,p4,dz1,dz2,dz real :: weight,hx,hy,hz,h,xd,yd,zd,xkern,r2,c(maxspec),ddx,ddy real :: rhoprof(2),rhoi real :: xl,yl,wx,wy,w real,parameter :: factor=.596831, hxmax=6.0, hymax=4.0, hzmax=150. ! For forward simulations, make a loop over the number of species; ! for backward simulations, make an additional loop over the ! releasepoints !*************************************************************************** do i=1,numpart if (itra1(i).ne.itime) goto 20 ! Determine age class of the particle itage=abs(itra1(i)-itramem(i)) do nage=1,nageclass if (itage.lt.lage(nage)) goto 33 end do 33 continue ! For special runs, interpolate the air density to the particle position !************************************************************************ !*********************************************************************** !AF IND_SOURCE switches between different units for concentrations at the source !Af NOTE that in backward simulations the release of particles takes place !Af at the receptor and the sampling at the source. !Af 1="mass" !Af 2="mass mixing ratio" !Af IND_RECEPTOR switches between different units for concentrations at the receptor !Af 1="mass" !Af 2="mass mixing ratio" !Af switches for the conccalcfile: !AF IND_SAMP = 0 : xmass * 1 !Af IND_SAMP = -1 : xmass / rho !Af ind_samp is defined in readcommand.f if ( ind_samp .eq. -1 ) then ix=int(xtra1(i)) jy=int(ytra1(i)) ixp=ix+1 jyp=jy+1 ddx=xtra1(i)-real(ix) ddy=ytra1(i)-real(jy) rddx=1.-ddx rddy=1.-ddy p1=rddx*rddy p2=ddx*rddy p3=rddx*ddy p4=ddx*ddy do il=2,nz if (height(il).gt.ztra1(i)) then indz=il-1 indzp=il goto 6 endif end do 6 continue dz1=ztra1(i)-height(indz) dz2=height(indzp)-ztra1(i) dz=1./(dz1+dz2) ! Take density from 2nd wind field in memory (accurate enough, no time interpolation needed) !***************************************************************************** do ind=indz,indzp rhoprof(ind-indz+1)=p1*rho(ix ,jy ,ind,2) & +p2*rho(ixp,jy ,ind,2) & +p3*rho(ix ,jyp,ind,2) & +p4*rho(ixp,jyp,ind,2) end do rhoi=(dz1*rhoprof(2)+dz2*rhoprof(1))*dz elseif (ind_samp.eq.0) then rhoi = 1. endif !**************************************************************************** ! 1. Evaluate grid concentrations using a uniform kernel of bandwidths dx, dy !**************************************************************************** ! For backward simulations, look from which release point the particle comes from ! For domain-filling trajectory option, npoint contains a consecutive particle ! number, not the release point information. Therefore, nrelpointer is set to 1 ! for the domain-filling option. !***************************************************************************** if ((ioutputforeachrelease.eq.0).or.(mdomainfill.eq.1)) then nrelpointer=1 else nrelpointer=npoint(i) endif do kz=1,numzgrid ! determine height of cell if (outheight(kz).gt.ztra1(i)) goto 21 end do 21 continue if (kz.le.numzgrid) then ! inside output domain !******************************** ! Do everything for mother domain !******************************** xl=(xtra1(i)*dx+xoutshift)/dxout yl=(ytra1(i)*dy+youtshift)/dyout ix=int(xl) if (xl.lt.0.) ix=ix-1 jy=int(yl) if (yl.lt.0.) jy=jy-1 ! if (i.eq.10000) write(*,*) itime,xtra1(i),ytra1(i),ztra1(i),xl,yl ! For particles aged less than 3 hours, attribute particle mass to grid cell ! it resides in rather than use the kernel, in order to avoid its smoothing effect. ! For older particles, use the uniform kernel. ! If a particle is close to the domain boundary, do not use the kernel either. !***************************************************************************** if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. & (xl.gt.real(numxgrid-1)-0.5).or. & (yl.gt.real(numygrid-1)-0.5)) then ! no kernel, direct attribution to grid cell if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgrid-1).and. & (jy.le.numygrid-1)) then do ks=1,nspec gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight end do endif else ! attribution via uniform kernel ddx=xl-real(ix) ! distance to left cell border ddy=yl-real(jy) ! distance to lower cell border if (ddx.gt.0.5) then ixp=ix+1 wx=1.5-ddx else ixp=ix-1 wx=0.5+ddx endif if (ddy.gt.0.5) then jyp=jy+1 wy=1.5-ddy else jyp=jy-1 wy=0.5+ddy endif ! Determine mass fractions for four grid points !********************************************** if ((ix.ge.0).and.(ix.le.numxgrid-1)) then if ((jy.ge.0).and.(jy.le.numygrid-1)) then w=wx*wy do ks=1,nspec gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & gridunc(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then w=wx*(1.-wy) do ks=1,nspec gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & gridunc(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif endif if ((ixp.ge.0).and.(ixp.le.numxgrid-1)) then if ((jyp.ge.0).and.(jyp.le.numygrid-1)) then w=(1.-wx)*(1.-wy) do ks=1,nspec gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & gridunc(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif if ((jy.ge.0).and.(jy.le.numygrid-1)) then w=(1.-wx)*wy do ks=1,nspec gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & gridunc(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif endif endif !************************************ ! Do everything for the nested domain !************************************ if (nested_output.eq.1) then xl=(xtra1(i)*dx+xoutshiftn)/dxoutn yl=(ytra1(i)*dy+youtshiftn)/dyoutn ix=int(xl) if (xl.lt.0.) ix=ix-1 jy=int(yl) if (yl.lt.0.) jy=jy-1 ! For particles aged less than 3 hours, attribute particle mass to grid cell ! it resides in rather than use the kernel, in order to avoid its smoothing effect. ! For older particles, use the uniform kernel. ! If a particle is close to the domain boundary, do not use the kernel either. !***************************************************************************** if ((itage.lt.10800).or.(xl.lt.0.5).or.(yl.lt.0.5).or. & (xl.gt.real(numxgridn-1)-0.5).or. & (yl.gt.real(numygridn-1)-0.5)) then ! no kernel, direct attribution to grid cell if ((ix.ge.0).and.(jy.ge.0).and.(ix.le.numxgridn-1).and. & (jy.le.numygridn-1)) then do ks=1,nspec griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight end do endif else ! attribution via uniform kernel ddx=xl-real(ix) ! distance to left cell border ddy=yl-real(jy) ! distance to lower cell border if (ddx.gt.0.5) then ixp=ix+1 wx=1.5-ddx else ixp=ix-1 wx=0.5+ddx endif if (ddy.gt.0.5) then jyp=jy+1 wy=1.5-ddy else jyp=jy-1 wy=0.5+ddy endif ! Determine mass fractions for four grid points !********************************************** if ((ix.ge.0).and.(ix.le.numxgridn-1)) then if ((jy.ge.0).and.(jy.le.numygridn-1)) then w=wx*wy do ks=1,nspec griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)= & griduncn(ix,jy,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then w=wx*(1.-wy) do ks=1,nspec griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)= & griduncn(ix,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif endif if ((ixp.ge.0).and.(ixp.le.numxgridn-1)) then if ((jyp.ge.0).and.(jyp.le.numygridn-1)) then w=(1.-wx)*(1.-wy) do ks=1,nspec griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)= & griduncn(ixp,jyp,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif if ((jy.ge.0).and.(jy.le.numygridn-1)) then w=(1.-wx)*wy do ks=1,nspec griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)= & griduncn(ixp,jy,kz,ks,nrelpointer,nclass(i),nage)+ & xmass1(i,ks)/rhoi*weight*w end do endif endif endif endif endif 20 continue end do !*********************************************************************** ! 2. Evaluate concentrations at receptor points, using the kernel method !*********************************************************************** do n=1,numreceptor ! Reset concentrations !********************* do ks=1,nspec c(ks)=0. end do ! Estimate concentration at receptor !*********************************** do i=1,numpart if (itra1(i).ne.itime) goto 40 itage=abs(itra1(i)-itramem(i)) hz=min(50.+0.3*sqrt(real(itage)),hzmax) zd=ztra1(i)/hz if (zd.gt.1.) goto 40 ! save computing time, leave loop hx=min((0.29+2.222e-3*sqrt(real(itage)))*dx+ & real(itage)*1.2e-5,hxmax) ! 80 km/day xd=(xtra1(i)-xreceptor(n))/hx if (xd*xd.gt.1.) goto 40 ! save computing time, leave loop hy=min((0.18+1.389e-3*sqrt(real(itage)))*dy+ & real(itage)*7.5e-6,hymax) ! 80 km/day yd=(ytra1(i)-yreceptor(n))/hy if (yd*yd.gt.1.) goto 40 ! save computing time, leave loop h=hx*hy*hz r2=xd*xd+yd*yd+zd*zd if (r2.lt.1.) then xkern=factor*(1.-r2) do ks=1,nspec c(ks)=c(ks)+xmass1(i,ks)*xkern/h end do endif 40 continue end do do ks=1,nspec creceptor(n,ks)=creceptor(n,ks)+2.*weight*c(ks)/receptorarea(n) end do end do end subroutine conccalc