source: trunk/src/boundcond_domainfill.f90 @ 28

Last change on this file since 28 was 4, checked in by mlanger, 11 years ago
File size: 20.7 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 boundcond_domainfill(itime,loutend)
23  !                                  i      i
24  !*****************************************************************************
25  !                                                                            *
26  ! Particles are created by this subroutine continuously throughout the       *
27  ! simulation at the boundaries of the domain-filling box.                    *
28  ! All particles carry the same amount of mass which alltogether comprises the*
29  ! mass of air within the box, which remains (more or less) constant.         *
30  !                                                                            *
31  !     Author: A. Stohl                                                       *
32  !                                                                            *
33  !     16 October 2002                                                        *
34  !                                                                            *
35  !*****************************************************************************
36  !                                                                            *
37  ! Variables:                                                                 *
38  !                                                                            *
39  ! nx_we(2)       grid indices for western and eastern boundary of domain-    *
40  !                filling trajectory calculations                             *
41  ! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
42  !                filling trajectory calculations                             *
43  !                                                                            *
44  !*****************************************************************************
45
46  use point_mod
47  use par_mod
48  use com_mod
49
50  implicit none
51
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
56
57  real :: windl(2),rhol(2)
58  real :: windhl(2),rhohl(2)
59  real :: windx,rhox
60  real :: deltaz,boundarea,fluxofmass
61
62  integer :: ixm,ixp,jym,jyp,indzm,mm
63  real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2),yh1(2)
64
65  integer :: idummy = -11
66
67
68  ! If domain-filling is global, no boundary conditions are needed
69  !***************************************************************
70
71  if (gdomainfill) return
72
73  accmasst=0.
74  numactiveparticles=0
75
76  ! Terminate trajectories that have left the domain, if domain-filling
77  ! trajectory calculation domain is not global
78  !********************************************************************
79
80  do i=1,numpart
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
87    endif
88    if (itra1(i).ne.-999999999) numactiveparticles= &
89         numactiveparticles+1
90  end do
91
92
93  ! Determine auxiliary variables for time interpolation
94  !*****************************************************
95
96  dt1=real(itime-memtime(1))
97  dt2=real(memtime(2)-itime)
98  dtt=1./(dt1+dt2)
99
100  ! Initialize auxiliary variable used to search for vacant storage space
101  !**********************************************************************
102
103  minpart=1
104
105  !***************************************
106  ! Western and eastern boundary condition
107  !***************************************
108
109  ! Loop from south to north
110  !*************************
111
112  do jy=ny_sn(1),ny_sn(2)
113
114  ! Loop over western (index 1) and eastern (index 2) boundary
115  !***********************************************************
116
117    do k=1,2
118
119  ! Loop over all release locations in a column
120  !********************************************
121
122      do j=1,numcolumn_we(k,jy)
123
124  ! Determine, for each release location, the area of the corresponding boundary
125  !*****************************************************************************
126
127        if (j.eq.1) then
128          deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2.
129        else if (j.eq.numcolumn_we(k,jy)) then
130  !        deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+
131  !    +        zcolumn_we(k,jy,j))/2.
132  ! In order to avoid taking a very high column for very many particles,
133  ! use the deltaz from one particle below instead
134          deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2.
135        else
136          deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2.
137        endif
138        if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then
139          boundarea=deltaz*111198.5/2.*dy
140        else
141          boundarea=deltaz*111198.5*dy
142        endif
143
144
145  ! Interpolate the wind velocity and density to the release location
146  !******************************************************************
147
148  ! Determine the model level below the release position
149  !*****************************************************
150
151        do i=2,nz
152          if (height(i).gt.zcolumn_we(k,jy,j)) then
153            indz=i-1
154            indzp=i
155            goto 6
156          endif
157        end do
1586       continue
159
160  ! Vertical distance to the level below and above current position
161  !****************************************************************
162
163        dz1=zcolumn_we(k,jy,j)-height(indz)
164        dz2=height(indzp)-zcolumn_we(k,jy,j)
165        dz=1./(dz1+dz2)
166
167  ! Vertical and temporal interpolation
168  !************************************
169
170        do m=1,2
171          indexh=memind(m)
172          do in=1,2
173            indzh=indz+in-1
174            windl(in)=uu(nx_we(k),jy,indzh,indexh)
175            rhol(in)=rho(nx_we(k),jy,indzh,indexh)
176          end do
177
178          windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
179          rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
180        end do
181
182        windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
183        rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
184
185  ! Calculate mass flux
186  !********************
187
188        fluxofmass=windx*rhox*boundarea*real(lsynctime)
189
190
191  ! If the mass flux is directed into the domain, add it to previous mass fluxes;
192  ! if it is out of the domain, set accumulated mass flux to zero
193  !******************************************************************************
194
195        if (k.eq.1) then
196          if (fluxofmass.ge.0.) then
197            acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass
198          else
199            acc_mass_we(k,jy,j)=0.
200          endif
201        else
202          if (fluxofmass.le.0.) then
203            acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass)
204          else
205            acc_mass_we(k,jy,j)=0.
206          endif
207        endif
208        accmasst=accmasst+acc_mass_we(k,jy,j)
209
210  ! If the accumulated mass exceeds half the mass that each particle shall carry,
211  ! one (or more) particle(s) is (are) released and the accumulated mass is
212  ! reduced by the mass of this (these) particle(s)
213  !******************************************************************************
214
215        if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then
216          mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ &
217               xmassperparticle)
218          acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- &
219               real(mmass)*xmassperparticle
220        else
221          mmass=0
222        endif
223
224        do m=1,mmass
225          do ipart=minpart,maxpart
226
227  ! If a vacant storage space is found, attribute everything to this array element
228  !*****************************************************************************
229
230            if (itra1(ipart).ne.itime) then
231
232  ! Assign particle positions
233  !**************************
234
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)
240              else
241                ytra1(ipart)=real(jy)+(ran1(idummy)-.5)
242              endif
243              if (j.eq.1) then
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.
249              else
250                ztra1(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* &
251                     (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))
252              endif
253
254  ! Interpolate PV to the particle position
255  !****************************************
256              ixm=int(xtra1(ipart))
257              jym=int(ytra1(ipart))
258              ixp=ixm+1
259              jyp=jym+1
260              ddx=xtra1(ipart)-real(ixm)
261              ddy=ytra1(ipart)-real(jym)
262              rddx=1.-ddx
263              rddy=1.-ddy
264              p1=rddx*rddy
265              p2=ddx*rddy
266              p3=rddx*ddy
267              p4=ddx*ddy
268              do i=2,nz
269                if (height(i).gt.ztra1(ipart)) then
270                  indzm=i-1
271                  indzp=i
272                  goto 26
273                endif
274              end do
27526            continue
276              dz1=ztra1(ipart)-height(indzm)
277              dz2=height(indzp)-ztra1(ipart)
278              dz=1./(dz1+dz2)
279              do mm=1,2
280                indexh=memind(mm)
281                do in=1,2
282                  indzh=indzm+in-1
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)
287                end do
288                yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
289              end do
290              pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
291              ylat=ylat0+ytra1(ipart)*dy
292              if (ylat.lt.0.) pvpart=-1.*pvpart
293
294
295  ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
296  !*****************************************************************************
297
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
305                idt(ipart)=mintime
306                itra1(ipart)=itime
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
312              else
313                goto 71
314              endif
315
316
317  ! Increase numpart, if necessary
318  !*******************************
319
320              numpart=max(numpart,ipart)
321              goto 73      ! Storage space has been found, stop searching
322            endif
323          end do
324          if (ipart.gt.maxpart) &
325               stop 'boundcond_domainfill.f: too many particles required'
32673        minpart=ipart+1
32771        continue
328        end do
329
330
331      end do
332    end do
333  end do
334
335
336  !*****************************************
337  ! Southern and northern boundary condition
338  !*****************************************
339
340  ! Loop from west to east
341  !***********************
342
343  do ix=nx_we(1),nx_we(2)
344
345  ! Loop over southern (index 1) and northern (index 2) boundary
346  !*************************************************************
347
348    do k=1,2
349      ylat=ylat0+real(ny_sn(k))*dy
350      cosfact=cos(ylat*pi180)
351
352  ! Loop over all release locations in a column
353  !********************************************
354
355      do j=1,numcolumn_sn(k,ix)
356
357  ! Determine, for each release location, the area of the corresponding boundary
358  !*****************************************************************************
359
360        if (j.eq.1) then
361          deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2.
362        else if (j.eq.numcolumn_sn(k,ix)) then
363  !        deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+
364  !    +        zcolumn_sn(k,ix,j))/2.
365  ! In order to avoid taking a very high column for very many particles,
366  ! use the deltaz from one particle below instead
367          deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2.
368        else
369          deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2.
370        endif
371        if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then
372          boundarea=deltaz*111198.5/2.*cosfact*dx
373        else
374          boundarea=deltaz*111198.5*cosfact*dx
375        endif
376
377
378  ! Interpolate the wind velocity and density to the release location
379  !******************************************************************
380
381  ! Determine the model level below the release position
382  !*****************************************************
383
384        do i=2,nz
385          if (height(i).gt.zcolumn_sn(k,ix,j)) then
386            indz=i-1
387            indzp=i
388            goto 16
389          endif
390        end do
39116      continue
392
393  ! Vertical distance to the level below and above current position
394  !****************************************************************
395
396        dz1=zcolumn_sn(k,ix,j)-height(indz)
397        dz2=height(indzp)-zcolumn_sn(k,ix,j)
398        dz=1./(dz1+dz2)
399
400  ! Vertical and temporal interpolation
401  !************************************
402
403        do m=1,2
404          indexh=memind(m)
405          do in=1,2
406            indzh=indz+in-1
407            windl(in)=vv(ix,ny_sn(k),indzh,indexh)
408            rhol(in)=rho(ix,ny_sn(k),indzh,indexh)
409          end do
410
411          windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
412          rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
413        end do
414
415        windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
416        rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
417
418  ! Calculate mass flux
419  !********************
420
421        fluxofmass=windx*rhox*boundarea*real(lsynctime)
422
423  ! If the mass flux is directed into the domain, add it to previous mass fluxes;
424  ! if it is out of the domain, set accumulated mass flux to zero
425  !******************************************************************************
426
427        if (k.eq.1) then
428          if (fluxofmass.ge.0.) then
429            acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass
430          else
431            acc_mass_sn(k,ix,j)=0.
432          endif
433        else
434          if (fluxofmass.le.0.) then
435            acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass)
436          else
437            acc_mass_sn(k,ix,j)=0.
438          endif
439        endif
440        accmasst=accmasst+acc_mass_sn(k,ix,j)
441
442  ! If the accumulated mass exceeds half the mass that each particle shall carry,
443  ! one (or more) particle(s) is (are) released and the accumulated mass is
444  ! reduced by the mass of this (these) particle(s)
445  !******************************************************************************
446
447        if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then
448          mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ &
449               xmassperparticle)
450          acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- &
451               real(mmass)*xmassperparticle
452        else
453          mmass=0
454        endif
455
456        do m=1,mmass
457          do ipart=minpart,maxpart
458
459  ! If a vacant storage space is found, attribute everything to this array element
460  !*****************************************************************************
461
462            if (itra1(ipart).ne.itime) then
463
464  ! Assign particle positions
465  !**************************
466
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)
472              else
473                xtra1(ipart)=real(ix)+(ran1(idummy)-.5)
474              endif
475              if (j.eq.1) then
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.
481              else
482                ztra1(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* &
483                     (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))
484              endif
485
486
487  ! Interpolate PV to the particle position
488  !****************************************
489              ixm=int(xtra1(ipart))
490              jym=int(ytra1(ipart))
491              ixp=ixm+1
492              jyp=jym+1
493              ddx=xtra1(ipart)-real(ixm)
494              ddy=ytra1(ipart)-real(jym)
495              rddx=1.-ddx
496              rddy=1.-ddy
497              p1=rddx*rddy
498              p2=ddx*rddy
499              p3=rddx*ddy
500              p4=ddx*ddy
501              do i=2,nz
502                if (height(i).gt.ztra1(ipart)) then
503                  indzm=i-1
504                  indzp=i
505                  goto 126
506                endif
507              end do
508126           continue
509              dz1=ztra1(ipart)-height(indzm)
510              dz2=height(indzp)-ztra1(ipart)
511              dz=1./(dz1+dz2)
512              do mm=1,2
513                indexh=memind(mm)
514                do in=1,2
515                  indzh=indzm+in-1
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)
520                end do
521                yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
522              end do
523              pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
524              if (ylat.lt.0.) pvpart=-1.*pvpart
525
526
527  ! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
528  !*****************************************************************************
529
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
537                idt(ipart)=mintime
538                itra1(ipart)=itime
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
544              else
545                goto 171
546              endif
547
548
549  ! Increase numpart, if necessary
550  !*******************************
551              numpart=max(numpart,ipart)
552              goto 173      ! Storage space has been found, stop searching
553            endif
554          end do
555          if (ipart.gt.maxpart) &
556               stop 'boundcond_domainfill.f: too many particles required'
557173       minpart=ipart+1
558171       continue
559        end do
560
561
562      end do
563    end do
564  end do
565
566
567  xm=0.
568  do i=1,numpart
569    if (itra1(i).eq.itime) xm=xm+xmass1(i,1)
570  end do
571
572  !write(*,*) itime,numactiveparticles,numparticlecount,numpart,
573  !    +xm,accmasst,xm+accmasst
574
575
576  ! If particles shall be dumped, then accumulated masses at the domain boundaries
577  ! must be dumped, too, to be used for later runs
578  !*****************************************************************************
579
580  if ((ipout.gt.0).and.(itime.eq.loutend)) then
581    open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
582         form='unformatted')
583    write(unitboundcond) numcolumn_we,numcolumn_sn, &
584         zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
585    close(unitboundcond)
586  endif
587
588end subroutine boundcond_domainfill
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG