source: flexpart.git/src/boundcond_domainfill_mpi.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 16b61a5, checked in by Espen Sollum ATMOS <eso@…>, 3 years ago

Reworked the domain-filling option (MPI). Fixed a slow loop which had errors in loop counter (MPI)

  • Property mode set to 100644
File size: 29.3 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!  CHANGES                                                                   
46!    08/2016 eso: MPI version:                                               
47!
48!   -Root process release particles and distributes to other processes.
49!    Temporary arrays are used, also for the non-root (receiving) processes.
50!   -The scheme can be improved by having all processes report numpart
51!    (keeping track of how many particles have left the domain), so that
52!    a proportional amount of new particles can be distributed (however
53!    we have a separate function called from timemanager that will
54!    redistribute particles among processes if there are imbalances)     
55!*****************************************************************************
56
57  use point_mod
58  use par_mod
59  use com_mod
60  use random_mod, only: ran1
61  use mpi_mod
62
63  implicit none
64
65  real :: dz,dz1,dz2,dt1,dt2,dtt,ylat,xm,cosfact,accmasst
66  integer :: itime,in,indz,indzp,i,loutend
67  integer :: j,k,ix,jy,m,indzh,indexh,minpart,ipart,mmass
68  integer :: numactiveparticles, numpart_total, rel_counter
69  integer,allocatable,dimension(:) ::  numrel_mpi !, numactiveparticles_mpi
70
71  real :: windl(2),rhol(2)
72  real :: windhl(2),rhohl(2)
73  real :: windx,rhox
74  real :: deltaz,boundarea,fluxofmass
75
76  integer :: ixm,ixp,jym,jyp,indzm,mm
77  real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2),yh1(2)
78
79  integer :: idummy = -11
80  integer :: mtag
81  logical :: first_call=.true.
82! Sizes of temporary arrays are maxpartfract*maxpart. Increase maxpartfract if
83! needed. 
84  real,parameter :: maxpartfract=0.1
85  integer :: tmp_size = int(maxpartfract*maxpart)
86
87! Use different seed for each process
88  if (first_call) then
89    idummy=idummy+mp_seed
90    first_call=.false.
91  end if
92
93
94! If domain-filling is global, no boundary conditions are needed
95!***************************************************************
96
97  if (gdomainfill) return
98
99  accmasst=0.
100  numactiveparticles=0
101! Keep track of active particles on each process
102  allocate(numrel_mpi(0:mp_partgroup_np-1))
103! numactiveparticles_mpi(0:mp_partgroup_np-1)
104
105! New particles to be released on each process
106  numrel_mpi(:)=0
107
108! Terminate trajectories that have left the domain, if domain-filling
109! trajectory calculation domain is not global. Done for all processes
110!********************************************************************
111
112  do i=1,numpart
113    if (itra1(i).eq.itime) then
114      if ((ytra1(i).gt.real(ny_sn(2))).or. &
115           (ytra1(i).lt.real(ny_sn(1)))) itra1(i)=-999999999
116      if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
117           ((xtra1(i).lt.real(nx_we(1))).or. &
118           (xtra1(i).gt.real(nx_we(2))))) itra1(i)=-999999999
119    endif
120    if (itra1(i).ne.-999999999) numactiveparticles= &
121         numactiveparticles+1
122  end do
123!  numactiveparticles_mpi(mp_partid) = numactiveparticles
124
125
126! Collect number of active particles from all processes
127  ! call MPI_Allgather(numactiveparticles, 1, MPI_INTEGER, &
128  !      &numactiveparticles_mpi, 1, MPI_INTEGER, mp_comm_used, mp_ierr)
129
130
131! Total number of new releases
132  numpart_total = 0
133
134
135! This section only done by root process
136!***************************************
137
138  if (lroot) then
139
140! Use separate arrays for newly released particles
141!*************************************************
142
143    allocate(itra1_tmp(tmp_size),npoint_tmp(tmp_size),nclass_tmp(tmp_size),&
144         & idt_tmp(tmp_size),itramem_tmp(tmp_size),itrasplit_tmp(tmp_size),&
145         & xtra1_tmp(tmp_size),ytra1_tmp(tmp_size),ztra1_tmp(tmp_size),&
146         & xmass1_tmp(tmp_size, maxspec))
147
148! Initialize all particles as non-existent   
149    itra1_tmp(:)=-999999999
150
151! Determine auxiliary variables for time interpolation
152!*****************************************************
153
154    dt1=real(itime-memtime(1))
155    dt2=real(memtime(2)-itime)
156    dtt=1./(dt1+dt2)
157
158! Initialize auxiliary variable used to search for vacant storage space
159!**********************************************************************
160
161    minpart=1
162
163!***************************************
164! Western and eastern boundary condition
165!***************************************
166
167! Loop from south to north
168!*************************
169
170    do jy=ny_sn(1),ny_sn(2)
171
172! Loop over western (index 1) and eastern (index 2) boundary
173!***********************************************************
174
175      do k=1,2
176
177! Loop over all release locations in a column
178!********************************************
179
180        do j=1,numcolumn_we(k,jy)
181
182! Determine, for each release location, the area of the corresponding boundary
183!*****************************************************************************
184
185          if (j.eq.1) then
186            deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2.
187          else if (j.eq.numcolumn_we(k,jy)) then
188!        deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+
189!    +        zcolumn_we(k,jy,j))/2.
190! In order to avoid taking a very high column for very many particles,
191! use the deltaz from one particle below instead
192            deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2.
193          else
194            deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2.
195          endif
196          if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then
197            boundarea=deltaz*111198.5/2.*dy
198          else
199            boundarea=deltaz*111198.5*dy
200          endif
201
202
203! Interpolate the wind velocity and density to the release location
204!******************************************************************
205
206! Determine the model level below the release position
207!*****************************************************
208
209          do i=2,nz
210            if (height(i).gt.zcolumn_we(k,jy,j)) then
211              indz=i-1
212              indzp=i
213              goto 6
214            endif
215          end do
2166         continue
217
218! Vertical distance to the level below and above current position
219!****************************************************************
220
221          dz1=zcolumn_we(k,jy,j)-height(indz)
222          dz2=height(indzp)-zcolumn_we(k,jy,j)
223          dz=1./(dz1+dz2)
224
225! Vertical and temporal interpolation
226!************************************
227
228          do m=1,2
229            indexh=memind(m)
230            do in=1,2
231              indzh=indz+in-1
232              windl(in)=uu(nx_we(k),jy,indzh,indexh)
233              rhol(in)=rho(nx_we(k),jy,indzh,indexh)
234            end do
235
236            windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
237            rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
238          end do
239
240          windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
241          rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
242
243! Calculate mass flux
244!********************
245
246          fluxofmass=windx*rhox*boundarea*real(lsynctime)
247
248
249! If the mass flux is directed into the domain, add it to previous mass fluxes;
250! if it is out of the domain, set accumulated mass flux to zero
251!******************************************************************************
252
253          if (k.eq.1) then
254            if (fluxofmass.ge.0.) then
255              acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass
256            else
257              acc_mass_we(k,jy,j)=0.
258            endif
259          else
260            if (fluxofmass.le.0.) then
261              acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass)
262            else
263              acc_mass_we(k,jy,j)=0.
264            endif
265          endif
266          accmasst=accmasst+acc_mass_we(k,jy,j)
267
268! If the accumulated mass exceeds half the mass that each particle shall carry,
269! one (or more) particle(s) is (are) released and the accumulated mass is
270! reduced by the mass of this (these) particle(s)
271!******************************************************************************
272
273          if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then
274            mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ &
275                 xmassperparticle)
276            acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- &
277                 real(mmass)*xmassperparticle
278          else
279            mmass=0
280          endif
281
282          do m=1,mmass
283            do ipart=minpart,maxpart
284
285! If a vacant storage space is found, attribute everything to this array element
286! TODO: for the MPI version this test can be removed, as all
287!       elements in _tmp arrays are initialized to zero
288!*****************************************************************************
289
290              if (itra1_tmp(ipart).ne.itime) then
291
292! Assign particle positions
293!**************************
294
295                xtra1_tmp(ipart)=real(nx_we(k))
296                if (jy.eq.ny_sn(1)) then
297                  ytra1_tmp(ipart)=real(jy)+0.5*ran1(idummy)
298                else if (jy.eq.ny_sn(2)) then
299                  ytra1_tmp(ipart)=real(jy)-0.5*ran1(idummy)
300                else
301                  ytra1_tmp(ipart)=real(jy)+(ran1(idummy)-.5)
302                endif
303                if (j.eq.1) then
304                  ztra1_tmp(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- &
305                       zcolumn_we(k,jy,1))/4.
306                else if (j.eq.numcolumn_we(k,jy)) then
307                  ztra1_tmp(ipart)=(2.*zcolumn_we(k,jy,j)+ &
308                       zcolumn_we(k,jy,j-1)+height(nz))/4.
309                else
310                  ztra1_tmp(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* &
311                       (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))
312                endif
313
314! Interpolate PV to the particle position
315!****************************************
316                ixm=int(xtra1_tmp(ipart))
317                jym=int(ytra1_tmp(ipart))
318                ixp=ixm+1
319                jyp=jym+1
320                ddx=xtra1_tmp(ipart)-real(ixm)
321                ddy=ytra1_tmp(ipart)-real(jym)
322                rddx=1.-ddx
323                rddy=1.-ddy
324                p1=rddx*rddy
325                p2=ddx*rddy
326                p3=rddx*ddy
327                p4=ddx*ddy
328                do i=2,nz
329                  if (height(i).gt.ztra1_tmp(ipart)) then
330                    indzm=i-1
331                    indzp=i
332                    goto 26
333                  endif
334                end do
33526              continue
336                dz1=ztra1_tmp(ipart)-height(indzm)
337                dz2=height(indzp)-ztra1_tmp(ipart)
338                dz=1./(dz1+dz2)
339                do mm=1,2
340                  indexh=memind(mm)
341                  do in=1,2
342                    indzh=indzm+in-1
343                    y1(in)=p1*pv(ixm,jym,indzh,indexh) &
344                         +p2*pv(ixp,jym,indzh,indexh) &
345                         +p3*pv(ixm,jyp,indzh,indexh) &
346                         +p4*pv(ixp,jyp,indzh,indexh)
347                  end do
348                  yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
349                end do
350                pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
351                ylat=ylat0+ytra1_tmp(ipart)*dy
352                if (ylat.lt.0.) pvpart=-1.*pvpart
353
354
355! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
356!*****************************************************************************
357
358                if (((ztra1_tmp(ipart).gt.3000.).and. &
359                     (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
360                  nclass_tmp(ipart)=min(int(ran1(idummy)* &
361                       real(nclassunc))+1,nclassunc)
362                  numactiveparticles=numactiveparticles+1
363                  numparticlecount=numparticlecount+1
364                  npoint_tmp(ipart)=numparticlecount
365                  idt_tmp(ipart)=mintime
366                  itra1_tmp(ipart)=itime
367                  itramem_tmp(ipart)=itra1_tmp(ipart)
368                  itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit
369                  xmass1_tmp(ipart,1)=xmassperparticle
370                  if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= &
371                       xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9
372                else
373                  goto 71
374                endif
375
376
377! Increase numpart, if necessary
378!*******************************
379
380                numpart_total=max(numpart_total,ipart)
381                goto 73      ! Storage space has been found, stop searching
382              endif
383            end do
384            if (ipart.gt.tmp_size) &
385                 stop 'boundcond_domainfill_mpi.f90: too many particles required'
38673          minpart=ipart+1
38771          continue
388          end do
389
390
391        end do
392      end do
393    end do
394
395
396!*****************************************
397! Southern and northern boundary condition
398!*****************************************
399
400! Loop from west to east
401!***********************
402
403    do ix=nx_we(1),nx_we(2)
404
405! Loop over southern (index 1) and northern (index 2) boundary
406!*************************************************************
407
408      do k=1,2
409        ylat=ylat0+real(ny_sn(k))*dy
410        cosfact=cos(ylat*pi180)
411
412! Loop over all release locations in a column
413!********************************************
414
415        do j=1,numcolumn_sn(k,ix)
416
417! Determine, for each release location, the area of the corresponding boundary
418!*****************************************************************************
419
420          if (j.eq.1) then
421            deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2.
422          else if (j.eq.numcolumn_sn(k,ix)) then
423!        deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+
424!    +        zcolumn_sn(k,ix,j))/2.
425! In order to avoid taking a very high column for very many particles,
426! use the deltaz from one particle below instead
427            deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2.
428          else
429            deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2.
430          endif
431          if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then
432            boundarea=deltaz*111198.5/2.*cosfact*dx
433          else
434            boundarea=deltaz*111198.5*cosfact*dx
435          endif
436
437
438! Interpolate the wind velocity and density to the release location
439!******************************************************************
440
441! Determine the model level below the release position
442!*****************************************************
443
444          do i=2,nz
445            if (height(i).gt.zcolumn_sn(k,ix,j)) then
446              indz=i-1
447              indzp=i
448              goto 16
449            endif
450          end do
45116        continue
452
453! Vertical distance to the level below and above current position
454!****************************************************************
455
456          dz1=zcolumn_sn(k,ix,j)-height(indz)
457          dz2=height(indzp)-zcolumn_sn(k,ix,j)
458          dz=1./(dz1+dz2)
459
460! Vertical and temporal interpolation
461!************************************
462
463          do m=1,2
464            indexh=memind(m)
465            do in=1,2
466              indzh=indz+in-1
467              windl(in)=vv(ix,ny_sn(k),indzh,indexh)
468              rhol(in)=rho(ix,ny_sn(k),indzh,indexh)
469            end do
470
471            windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
472            rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
473          end do
474
475          windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
476          rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
477
478! Calculate mass flux
479!********************
480
481          fluxofmass=windx*rhox*boundarea*real(lsynctime)
482
483! If the mass flux is directed into the domain, add it to previous mass fluxes;
484! if it is out of the domain, set accumulated mass flux to zero
485!******************************************************************************
486
487          if (k.eq.1) then
488            if (fluxofmass.ge.0.) then
489              acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass
490            else
491              acc_mass_sn(k,ix,j)=0.
492            endif
493          else
494            if (fluxofmass.le.0.) then
495              acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass)
496            else
497              acc_mass_sn(k,ix,j)=0.
498            endif
499          endif
500          accmasst=accmasst+acc_mass_sn(k,ix,j)
501
502! If the accumulated mass exceeds half the mass that each particle shall carry,
503! one (or more) particle(s) is (are) released and the accumulated mass is
504! reduced by the mass of this (these) particle(s)
505!******************************************************************************
506
507          if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then
508            mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ &
509                 xmassperparticle)
510            acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- &
511                 real(mmass)*xmassperparticle
512          else
513            mmass=0
514          endif
515
516          do m=1,mmass
517            do ipart=minpart,maxpart
518
519! If a vacant storage space is found, attribute everything to this array element
520!*****************************************************************************
521
522              if (itra1_tmp(ipart).ne.itime) then
523
524! Assign particle positions
525!**************************
526
527                ytra1_tmp(ipart)=real(ny_sn(k))
528                if (ix.eq.nx_we(1)) then
529                  xtra1_tmp(ipart)=real(ix)+0.5*ran1(idummy)
530                else if (ix.eq.nx_we(2)) then
531                  xtra1_tmp(ipart)=real(ix)-0.5*ran1(idummy)
532                else
533                  xtra1_tmp(ipart)=real(ix)+(ran1(idummy)-.5)
534                endif
535                if (j.eq.1) then
536                  ztra1_tmp(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- &
537                       zcolumn_sn(k,ix,1))/4.
538                else if (j.eq.numcolumn_sn(k,ix)) then
539                  ztra1_tmp(ipart)=(2.*zcolumn_sn(k,ix,j)+ &
540                       zcolumn_sn(k,ix,j-1)+height(nz))/4.
541                else
542                  ztra1_tmp(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* &
543                       (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))
544                endif
545
546
547! Interpolate PV to the particle position
548!****************************************
549                ixm=int(xtra1_tmp(ipart))
550                jym=int(ytra1_tmp(ipart))
551                ixp=ixm+1
552                jyp=jym+1
553                ddx=xtra1_tmp(ipart)-real(ixm)
554                ddy=ytra1_tmp(ipart)-real(jym)
555                rddx=1.-ddx
556                rddy=1.-ddy
557                p1=rddx*rddy
558                p2=ddx*rddy
559                p3=rddx*ddy
560                p4=ddx*ddy
561                do i=2,nz
562                  if (height(i).gt.ztra1_tmp(ipart)) then
563                    indzm=i-1
564                    indzp=i
565                    goto 126
566                  endif
567                end do
568126             continue
569                dz1=ztra1_tmp(ipart)-height(indzm)
570                dz2=height(indzp)-ztra1_tmp(ipart)
571                dz=1./(dz1+dz2)
572                do mm=1,2
573                  indexh=memind(mm)
574                  do in=1,2
575                    indzh=indzm+in-1
576                    y1(in)=p1*pv(ixm,jym,indzh,indexh) &
577                         +p2*pv(ixp,jym,indzh,indexh) &
578                         +p3*pv(ixm,jyp,indzh,indexh) &
579                         +p4*pv(ixp,jyp,indzh,indexh)
580                  end do
581                  yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
582                end do
583                pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
584                if (ylat.lt.0.) pvpart=-1.*pvpart
585
586
587! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
588!*****************************************************************************
589
590                if (((ztra1_tmp(ipart).gt.3000.).and. &
591                     (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
592                  nclass_tmp(ipart)=min(int(ran1(idummy)* &
593                       real(nclassunc))+1,nclassunc)
594                  numactiveparticles=numactiveparticles+1
595                  numparticlecount=numparticlecount+1
596                  npoint_tmp(ipart)=numparticlecount
597                  idt_tmp(ipart)=mintime
598                  itra1_tmp(ipart)=itime
599                  itramem_tmp(ipart)=itra1_tmp(ipart)
600                  itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit
601                  xmass1_tmp(ipart,1)=xmassperparticle
602                  if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= &
603                       xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9
604                else
605                  goto 171
606                endif
607
608
609! Increase numpart, if necessary
610!*******************************
611                numpart_total=max(numpart_total,ipart)
612                goto 173      ! Storage space has been found, stop searching
613              endif
614            end do
615            if (ipart.gt.tmp_size) &
616                 stop 'boundcond_domainfill.f: too many particles required'
617173         minpart=ipart+1
618171         continue
619          end do
620
621
622        end do
623      end do
624    end do
625
626
627! xm=0.
628! do i=1,numpart_total
629!   if (itra1_tmp(i).eq.itime) xm=xm+xmass1(i,1)
630! end do
631
632!write(*,*) itime,numactiveparticles,numparticlecount,numpart,
633!    +xm,accmasst,xm+accmasst
634
635  end if ! if lroot
636
637! Distribute the number of particles to be released
638! *************************************************
639  call MPI_Bcast(numpart_total, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)
640
641  do i=0, mp_partgroup_np-1
642    numrel_mpi(i) = numpart_total/mp_partgroup_np
643    if (i.lt.mod(numpart_total,mp_partgroup_np)) numrel_mpi(i) = numrel_mpi(i) + 1
644  end do
645
646! Allocate temporary arrays for receiving processes
647  if (.not.lroot) then
648    allocate(itra1_tmp(numrel_mpi(mp_partid)),&
649         & npoint_tmp(numrel_mpi(mp_partid)),&
650         & nclass_tmp(numrel_mpi(mp_partid)),&
651         & idt_tmp(numrel_mpi(mp_partid)),&
652         & itramem_tmp(numrel_mpi(mp_partid)),&
653         & itrasplit_tmp(numrel_mpi(mp_partid)),&
654         & xtra1_tmp(numrel_mpi(mp_partid)),&
655         & ytra1_tmp(numrel_mpi(mp_partid)),&
656         & ztra1_tmp(numrel_mpi(mp_partid)),&
657         & xmass1_tmp(numrel_mpi(mp_partid),maxspec))
658     
659! Initialize all particles as non-existent   
660    itra1_tmp(:)=-999999999
661  end if
662
663! Distribute particles
664! Keep track of released particles so far
665  rel_counter = 0
666  mtag = 1000
667
668  do i=0, mp_partgroup_np-1
669
670! For root process, nothing to do except update release count
671    if (i.eq.0) then
672      rel_counter = rel_counter + numrel_mpi(i)
673      cycle
674    end if
675     
676! Send particles from root to non-root processes
677    if (lroot.and.numrel_mpi(i).gt.0) then
678
679      call MPI_SEND(nclass_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
680           &numrel_mpi(i),MPI_INTEGER,i,mtag+1*i,mp_comm_used,mp_ierr)
681
682      call MPI_SEND(npoint_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
683           &numrel_mpi(i),MPI_INTEGER,i,mtag+2*i,mp_comm_used,mp_ierr)
684
685      call MPI_SEND(itra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
686           &numrel_mpi(i),MPI_INTEGER,i,mtag+3*i,mp_comm_used,mp_ierr)
687
688      call MPI_SEND(idt_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
689           &numrel_mpi(i),MPI_INTEGER,i,mtag+4*i,mp_comm_used,mp_ierr)
690
691      call MPI_SEND(itramem_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
692           &numrel_mpi(i),MPI_INTEGER,i,mtag+5*i,mp_comm_used,mp_ierr)
693
694      call MPI_SEND(itrasplit_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
695           &numrel_mpi(i),MPI_INTEGER,i,mtag+6*i,mp_comm_used,mp_ierr)
696
697      call MPI_SEND(xtra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
698           &numrel_mpi(i),mp_dp,i,mtag+7*i,mp_comm_used,mp_ierr)
699
700      call MPI_SEND(ytra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
701           &numrel_mpi(i),mp_dp,i,mtag+8*i,mp_comm_used,mp_ierr)
702
703      call MPI_SEND(ztra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
704           &numrel_mpi(i),mp_sp,i,mtag+9*i,mp_comm_used,mp_ierr)
705
706      do j=1,nspec
707        call MPI_SEND(xmass1_tmp(rel_counter+1:rel_counter+numrel_mpi(i),j),&
708             &numrel_mpi(i),mp_sp,i,mtag+(9+j)*i,mp_comm_used,mp_ierr)
709      end do
710
711! Non-root processes issue receive requests
712    else if (i.eq.mp_partid.and.numrel_mpi(i).gt.0) then
713      call MPI_RECV(nclass_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
714           &MPI_INTEGER,id_root,mtag+1*i,mp_comm_used,mp_status,mp_ierr)
715
716      call MPI_RECV(npoint_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
717           &MPI_INTEGER,id_root,mtag+2*i,mp_comm_used,mp_status,mp_ierr)
718
719      call MPI_RECV(itra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
720           &MPI_INTEGER,id_root,mtag+3*i,mp_comm_used,mp_status,mp_ierr)
721
722      call MPI_RECV(idt_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
723           &MPI_INTEGER,id_root,mtag+4*i,mp_comm_used,mp_status,mp_ierr)
724
725      call MPI_RECV(itramem_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
726           &MPI_INTEGER,id_root,mtag+5*i,mp_comm_used,mp_status,mp_ierr)
727
728      call MPI_RECV(itrasplit_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
729           &MPI_INTEGER,id_root,mtag+6*i,mp_comm_used,mp_status,mp_ierr)
730
731      call MPI_RECV(xtra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
732           &mp_dp,id_root,mtag+7*i,mp_comm_used,mp_status,mp_ierr)
733
734      call MPI_RECV(ytra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
735           &mp_dp,id_root,mtag+8*i,mp_comm_used,mp_status,mp_ierr)
736
737      call MPI_RECV(ztra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
738           &mp_sp,id_root,mtag+9*i,mp_comm_used,mp_status,mp_ierr)
739
740      do j=1,nspec
741        call MPI_RECV(xmass1_tmp(1:numrel_mpi(i),j),numrel_mpi(i),&
742             &mp_sp,id_root,mtag+(9+j)*i,mp_comm_used,mp_status,mp_ierr)
743
744      end do
745    end if
746    rel_counter = rel_counter + numrel_mpi(i)
747  end do
748
749! Find free storage space for the new particles.
750! This section is independent of the redistribution scheme used
751! ********************************************************************
752
753! Keep track of released particles so far
754  minpart=1
755   
756! The algorithm should be correct also for root process
757  do i=1, numrel_mpi(mp_partid)
758    do ipart=minpart, maxpart
759      if (itra1(ipart).ne.itime) then
760        itra1(ipart) = itra1_tmp(i)
761        npoint(ipart) = npoint_tmp(i)
762        nclass(ipart) = nclass_tmp(i)
763        idt(ipart) = idt_tmp(i)
764        itramem(ipart) = itramem_tmp(i)
765        itrasplit(ipart) = itrasplit_tmp(i)
766        xtra1(ipart) = xtra1_tmp(i)
767        ytra1(ipart) = ytra1_tmp(i)
768        ztra1(ipart) = ztra1_tmp(i)
769        xmass1(ipart,:) = xmass1_tmp(i,:)
770! Increase numpart, if necessary
771        numpart=max(numpart,ipart)
772        goto 200 ! Storage space has been found, stop searching
773      end if
774    end do
775200 minpart=ipart+1
776  end do
777
778! If particles shall be dumped, then accumulated masses at the domain boundaries
779! must be dumped, too, to be used for later runs
780!*****************************************************************************
781
782  if ((ipout.gt.0).and.(itime.eq.loutend)) then
783    if (lroot) then
784      call mpif_mtime('iotime',0)
785      open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
786           form='unformatted')
787      write(unitboundcond) numcolumn_we,numcolumn_sn, &
788           zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
789      close(unitboundcond)
790      call mpif_mtime('iotime',1)
791    end if
792  endif
793
794! Deallocate temporary arrays
795  deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,&
796       & xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp,numrel_mpi)
797! numactiveparticles_mpi
798
799
800end subroutine boundcond_domainfill
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG