source: flexpart.git/src/boundcond_domainfill.f90 @ 38b7917

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since 38b7917 was 8a65cb0, checked in by Espen Sollum ATMOS <espen@…>, 9 years ago

Added code, makefile for dev branch

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