source: flexpart.git/src/boundcond_domainfill_mpi.f90 @ d404d98

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie
Last change on this file since d404d98 was 0f7835d, checked in by Espen Sollum ATMOS <eso@…>, 8 years ago

Changed a variable in domain-filling option (MPI). Deleted an unused function from mpi_mod.f90

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