52 real :: dz,dz1,dz2,
ran1,dt1,dt2,dtt,ylat,xm,cosfact,accmasst
53 integer :: itime,in,indz,indzp,i,loutend
54 integer :: j,k,ix,jy,m,indzh,indexh,minpart,ipart,mmass
55 integer :: numactiveparticles
57 real :: windl(2),rhol(2)
58 real :: windhl(2),rhohl(2)
60 real :: deltaz,boundarea,fluxofmass
62 integer :: ixm,ixp,jym,jyp,indzm,mm
63 real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2),yh1(2)
65 integer :: idummy = -11
71 if (gdomainfill)
return
81 if (itra1(i).eq.itime)
then
82 if ((ytra1(i).gt.
real(ny_sn(2))).or. &
83 (ytra1(i).lt.
real(ny_sn(1)))) itra1(i)=-999999999
84 if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
85 ((xtra1(i).lt.
real(nx_we(1))).or. &
86 (xtra1(i).gt.
real(nx_we(2))))) itra1(i)=-999999999
88 if (itra1(i).ne.-999999999) numactiveparticles= &
96 dt1=
real(itime-memtime(1))
97 dt2=
real(memtime(2)-itime)
112 do jy=ny_sn(1),ny_sn(2)
122 do j=1,numcolumn_we(k,jy)
128 deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2.
129 else if (j.eq.numcolumn_we(k,jy))
then
134 deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2.
136 deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2.
138 if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2)))
then
139 boundarea=deltaz*111198.5/2.*dy
141 boundarea=deltaz*111198.5*dy
152 if (height(i).gt.zcolumn_we(k,jy,j))
then
163 dz1=zcolumn_we(k,jy,j)-height(indz)
164 dz2=height(indzp)-zcolumn_we(k,jy,j)
174 windl(in)=uu(nx_we(k),jy,indzh,indexh)
175 rhol(in)=rho(nx_we(k),jy,indzh,indexh)
178 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
179 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
182 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
183 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
188 fluxofmass=windx*rhox*boundarea*
real(lsynctime)
196 if (fluxofmass.ge.0.)
then
197 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass
199 acc_mass_we(k,jy,j)=0.
202 if (fluxofmass.le.0.)
then
203 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass)
205 acc_mass_we(k,jy,j)=0.
208 accmasst=accmasst+acc_mass_we(k,jy,j)
215 if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.)
then
216 mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ &
218 acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- &
219 real(mmass)*xmassperparticle
225 do ipart=minpart,maxpart
230 if (itra1(ipart).ne.itime)
then
235 xtra1(ipart)=
real(nx_we(k))
236 if (jy.eq.ny_sn(1))
then
237 ytra1(ipart)=
real(jy)+0.5*
ran1(idummy)
238 else if (jy.eq.ny_sn(2))
then
239 ytra1(ipart)=
real(jy)-0.5*
ran1(idummy)
241 ytra1(ipart)=
real(jy)+(
ran1(idummy)-.5)
244 ztra1(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- &
245 zcolumn_we(k,jy,1))/4.
246 else if (j.eq.numcolumn_we(k,jy))
then
247 ztra1(ipart)=(2.*zcolumn_we(k,jy,j)+ &
248 zcolumn_we(k,jy,j-1)+height(nz))/4.
250 ztra1(ipart)=zcolumn_we(k,jy,j-1)+
ran1(idummy)* &
251 (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))
256 ixm=int(xtra1(ipart))
257 jym=int(ytra1(ipart))
260 ddx=xtra1(ipart)-
real(ixm)
261 ddy=ytra1(ipart)-
real(jym)
269 if (height(i).gt.ztra1(ipart))
then
276 dz1=ztra1(ipart)-height(indzm)
277 dz2=height(indzp)-ztra1(ipart)
283 y1(in)=p1*pv(ixm,jym,indzh,indexh) &
284 +p2*pv(ixp,jym,indzh,indexh) &
285 +p3*pv(ixm,jyp,indzh,indexh) &
286 +p4*pv(ixp,jyp,indzh,indexh)
288 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
290 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
291 ylat=ylat0+ytra1(ipart)*dy
292 if (ylat.lt.0.) pvpart=-1.*pvpart
298 if (((ztra1(ipart).gt.3000.).and. &
299 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1))
then
300 nclass(ipart)=min(int(
ran1(idummy)* &
301 real(nclassunc))+1,nclassunc)
302 numactiveparticles=numactiveparticles+1
303 numparticlecount=numparticlecount+1
304 npoint(ipart)=numparticlecount
307 itramem(ipart)=itra1(ipart)
308 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
309 xmass1(ipart,1)=xmassperparticle
310 if (mdomainfill.eq.2) xmass1(ipart,1)= &
311 xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9
320 numpart=max(numpart,ipart)
324 if (ipart.gt.maxpart) &
325 stop
'boundcond_domainfill.f: too many particles required'
343 do ix=nx_we(1),nx_we(2)
349 ylat=ylat0+
real(ny_sn(k))*dy
350 cosfact=cos(ylat*pi180)
355 do j=1,numcolumn_sn(k,ix)
361 deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2.
362 else if (j.eq.numcolumn_sn(k,ix))
then
367 deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2.
369 deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2.
371 if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2)))
then
372 boundarea=deltaz*111198.5/2.*cosfact*dx
374 boundarea=deltaz*111198.5*cosfact*dx
385 if (height(i).gt.zcolumn_sn(k,ix,j))
then
396 dz1=zcolumn_sn(k,ix,j)-height(indz)
397 dz2=height(indzp)-zcolumn_sn(k,ix,j)
407 windl(in)=vv(ix,ny_sn(k),indzh,indexh)
408 rhol(in)=rho(ix,ny_sn(k),indzh,indexh)
411 windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
412 rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
415 windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
416 rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
421 fluxofmass=windx*rhox*boundarea*
real(lsynctime)
428 if (fluxofmass.ge.0.)
then
429 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass
431 acc_mass_sn(k,ix,j)=0.
434 if (fluxofmass.le.0.)
then
435 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass)
437 acc_mass_sn(k,ix,j)=0.
440 accmasst=accmasst+acc_mass_sn(k,ix,j)
447 if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.)
then
448 mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ &
450 acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- &
451 real(mmass)*xmassperparticle
457 do ipart=minpart,maxpart
462 if (itra1(ipart).ne.itime)
then
467 ytra1(ipart)=
real(ny_sn(k))
468 if (ix.eq.nx_we(1))
then
469 xtra1(ipart)=
real(ix)+0.5*
ran1(idummy)
470 else if (ix.eq.nx_we(2))
then
471 xtra1(ipart)=
real(ix)-0.5*
ran1(idummy)
473 xtra1(ipart)=
real(ix)+(
ran1(idummy)-.5)
476 ztra1(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- &
477 zcolumn_sn(k,ix,1))/4.
478 else if (j.eq.numcolumn_sn(k,ix))
then
479 ztra1(ipart)=(2.*zcolumn_sn(k,ix,j)+ &
480 zcolumn_sn(k,ix,j-1)+height(nz))/4.
482 ztra1(ipart)=zcolumn_sn(k,ix,j-1)+
ran1(idummy)* &
483 (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))
489 ixm=int(xtra1(ipart))
490 jym=int(ytra1(ipart))
493 ddx=xtra1(ipart)-
real(ixm)
494 ddy=ytra1(ipart)-
real(jym)
502 if (height(i).gt.ztra1(ipart))
then
509 dz1=ztra1(ipart)-height(indzm)
510 dz2=height(indzp)-ztra1(ipart)
516 y1(in)=p1*pv(ixm,jym,indzh,indexh) &
517 +p2*pv(ixp,jym,indzh,indexh) &
518 +p3*pv(ixm,jyp,indzh,indexh) &
519 +p4*pv(ixp,jyp,indzh,indexh)
521 yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
523 pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
524 if (ylat.lt.0.) pvpart=-1.*pvpart
530 if (((ztra1(ipart).gt.3000.).and. &
531 (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1))
then
532 nclass(ipart)=min(int(
ran1(idummy)* &
533 real(nclassunc))+1,nclassunc)
534 numactiveparticles=numactiveparticles+1
535 numparticlecount=numparticlecount+1
536 npoint(ipart)=numparticlecount
539 itramem(ipart)=itra1(ipart)
540 itrasplit(ipart)=itra1(ipart)+ldirect*itsplit
541 xmass1(ipart,1)=xmassperparticle
542 if (mdomainfill.eq.2) xmass1(ipart,1)= &
543 xmass1(ipart,1)*pvpart*48./29.*ozonescale/10.**9
551 numpart=max(numpart,ipart)
555 if (ipart.gt.maxpart) &
556 stop
'boundcond_domainfill.f: too many particles required'
569 if (itra1(i).eq.itime) xm=xm+xmass1(i,1)
580 if ((ipout.gt.0).and.(itime.eq.loutend))
then
581 open(unitboundcond,file=path(2)(1:length(2))//
'boundcond.bin', &
583 write(unitboundcond) numcolumn_we,numcolumn_sn, &
584 zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
subroutine boundcond_domainfill(itime, loutend)