boundcond_domainfill_mpi.f90 28 KB
Newer Older
1
2
! SPDX-FileCopyrightText: FLEXPART 1998-2019, see flexpart_license.txt
! SPDX-License-Identifier: GPL-3.0-or-later
3

4
subroutine boundcond_domainfill(itime,loutend)
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
!                                  i      i
!*****************************************************************************
!                                                                            *
! Particles are created by this subroutine continuously throughout the       *
! simulation at the boundaries of the domain-filling box.                    *
! All particles carry the same amount of mass which alltogether comprises the*
! mass of air within the box, which remains (more or less) constant.         *
!                                                                            *
!     Author: A. Stohl                                                       *
!                                                                            *
!     16 October 2002                                                        *
!                                                                            *
!*****************************************************************************
!                                                                            *
! Variables:                                                                 *
!                                                                            *
! nx_we(2)       grid indices for western and eastern boundary of domain-    *
!                filling trajectory calculations                             *
! ny_sn(2)       grid indices for southern and northern boundary of domain-  *
!                filling trajectory calculations                             *
!                                                                            *
!*****************************************************************************
!  CHANGES                                                                   
!    08/2016 eso: MPI version:                                                
!
!   -Root process release particles and distributes to other processes.
!    Temporary arrays are used, also for the non-root (receiving) processes.
!   -The scheme can be improved by having all processes report numpart
!    (keeping track of how many particles have left the domain), so that
!    a proportional amount of new particles can be distributed (however
!    we have a separate function called from timemanager that will
!    redistribute particles among processes if there are imbalances)     
!*****************************************************************************
38
39
40
41
42
43
44
45
46
47
48
49

  use point_mod
  use par_mod
  use com_mod
  use random_mod, only: ran1
  use mpi_mod

  implicit none

  real :: dz,dz1,dz2,dt1,dt2,dtt,ylat,xm,cosfact,accmasst
  integer :: itime,in,indz,indzp,i,loutend
  integer :: j,k,ix,jy,m,indzh,indexh,minpart,ipart,mmass
50
51
  integer :: numactiveparticles, numpart_total, rel_counter
  integer,allocatable,dimension(:) ::  numrel_mpi !, numactiveparticles_mpi
52
53
54
55
56
57
58
59
60
61

  real :: windl(2),rhol(2)
  real :: windhl(2),rhohl(2)
  real :: windx,rhox
  real :: deltaz,boundarea,fluxofmass

  integer :: ixm,ixp,jym,jyp,indzm,mm
  real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2),yh1(2)

  integer :: idummy = -11
62
  integer :: mtag
63
  logical :: first_call=.true.
64
65
66
67
! Sizes of temporary arrays are maxpartfract*maxpart. Increase maxpartfract if 
! needed.  
  real,parameter :: maxpartfract=0.1
  integer :: tmp_size = int(maxpartfract*maxpart)
68

69
! Use different seed for each process
70
71
72
73
74
75
  if (first_call) then
    idummy=idummy+mp_seed
    first_call=.false.
  end if


76
77
! If domain-filling is global, no boundary conditions are needed
!***************************************************************
78
79
80
81
82

  if (gdomainfill) return

  accmasst=0.
  numactiveparticles=0
83
84
85
! Keep track of active particles on each process
  allocate(numrel_mpi(0:mp_partgroup_np-1))
! numactiveparticles_mpi(0:mp_partgroup_np-1)
86

87
88
89
90
91
92
! New particles to be released on each process
  numrel_mpi(:)=0

! Terminate trajectories that have left the domain, if domain-filling
! trajectory calculation domain is not global. Done for all processes
!********************************************************************
93
94
95
96
97
98
99
100
101
102
103
104

  do i=1,numpart
    if (itra1(i).eq.itime) then
      if ((ytra1(i).gt.real(ny_sn(2))).or. &
           (ytra1(i).lt.real(ny_sn(1)))) itra1(i)=-999999999
      if (((.not.xglobal).or.(nx_we(2).ne.(nx-2))).and. &
           ((xtra1(i).lt.real(nx_we(1))).or. &
           (xtra1(i).gt.real(nx_we(2))))) itra1(i)=-999999999
    endif
    if (itra1(i).ne.-999999999) numactiveparticles= &
         numactiveparticles+1
  end do
105
!  numactiveparticles_mpi(mp_partid) = numactiveparticles
106
107


108
109
110
! Collect number of active particles from all processes
  ! call MPI_Allgather(numactiveparticles, 1, MPI_INTEGER, &
  !      &numactiveparticles_mpi, 1, MPI_INTEGER, mp_comm_used, mp_ierr)
111
112


113
114
! Total number of new releases
  numpart_total = 0
115
116


117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
! This section only done by root process
!***************************************

  if (lroot) then

! Use separate arrays for newly released particles
!*************************************************

    allocate(itra1_tmp(tmp_size),npoint_tmp(tmp_size),nclass_tmp(tmp_size),&
         & idt_tmp(tmp_size),itramem_tmp(tmp_size),itrasplit_tmp(tmp_size),&
         & xtra1_tmp(tmp_size),ytra1_tmp(tmp_size),ztra1_tmp(tmp_size),&
         & xmass1_tmp(tmp_size, maxspec))

! Initialize all particles as non-existent    
    itra1_tmp(:)=-999999999

! Determine auxiliary variables for time interpolation
!*****************************************************
135

136
137
138
    dt1=real(itime-memtime(1))
    dt2=real(memtime(2)-itime)
    dtt=1./(dt1+dt2)
139

140
141
! Initialize auxiliary variable used to search for vacant storage space
!**********************************************************************
142

143
    minpart=1
144

145
146
147
!***************************************
! Western and eastern boundary condition
!***************************************
148

149
150
! Loop from south to north
!*************************
151

152
    do jy=ny_sn(1),ny_sn(2)
153

154
155
! Loop over western (index 1) and eastern (index 2) boundary
!***********************************************************
156

157
      do k=1,2
158

159
160
! Loop over all release locations in a column
!********************************************
161

162
        do j=1,numcolumn_we(k,jy)
163

164
165
! Determine, for each release location, the area of the corresponding boundary
!*****************************************************************************
166

167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
          if (j.eq.1) then
            deltaz=(zcolumn_we(k,jy,2)+zcolumn_we(k,jy,1))/2.
          else if (j.eq.numcolumn_we(k,jy)) then
!        deltaz=height(nz)-(zcolumn_we(k,jy,j-1)+
!    +        zcolumn_we(k,jy,j))/2.
! In order to avoid taking a very high column for very many particles,
! use the deltaz from one particle below instead
            deltaz=(zcolumn_we(k,jy,j)-zcolumn_we(k,jy,j-2))/2.
          else
            deltaz=(zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))/2.
          endif
          if ((jy.eq.ny_sn(1)).or.(jy.eq.ny_sn(2))) then
            boundarea=deltaz*111198.5/2.*dy
          else
            boundarea=deltaz*111198.5*dy
182
183
184
          endif


185
186
! Interpolate the wind velocity and density to the release location
!******************************************************************
187

188
189
! Determine the model level below the release position
!*****************************************************
190

191
192
193
194
195
196
          do i=2,nz
            if (height(i).gt.zcolumn_we(k,jy,j)) then
              indz=i-1
              indzp=i
              goto 6
            endif
197
          end do
198
6         continue
199

200
201
! Vertical distance to the level below and above current position
!****************************************************************
202

203
204
205
          dz1=zcolumn_we(k,jy,j)-height(indz)
          dz2=height(indzp)-zcolumn_we(k,jy,j)
          dz=1./(dz1+dz2)
206

207
208
! Vertical and temporal interpolation
!************************************
209

210
211
212
213
214
215
216
          do m=1,2
            indexh=memind(m)
            do in=1,2
              indzh=indz+in-1
              windl(in)=uu(nx_we(k),jy,indzh,indexh)
              rhol(in)=rho(nx_we(k),jy,indzh,indexh)
            end do
217

218
219
220
221
222
223
224
225
226
227
228
            windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
            rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
          end do

          windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
          rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt

! Calculate mass flux
!********************

          fluxofmass=windx*rhox*boundarea*real(lsynctime)
229
230


231
232
233
234
235
236
237
238
239
240
! If the mass flux is directed into the domain, add it to previous mass fluxes;
! if it is out of the domain, set accumulated mass flux to zero
!******************************************************************************

          if (k.eq.1) then
            if (fluxofmass.ge.0.) then
              acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+fluxofmass
            else
              acc_mass_we(k,jy,j)=0.
            endif
241
          else
242
243
244
245
246
            if (fluxofmass.le.0.) then
              acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)+abs(fluxofmass)
            else
              acc_mass_we(k,jy,j)=0.
            endif
247
          endif
248
249
250
251
252
253
254
255
256
257
258
259
          accmasst=accmasst+acc_mass_we(k,jy,j)

! If the accumulated mass exceeds half the mass that each particle shall carry,
! one (or more) particle(s) is (are) released and the accumulated mass is
! reduced by the mass of this (these) particle(s)
!******************************************************************************

          if (acc_mass_we(k,jy,j).ge.xmassperparticle/2.) then
            mmass=int((acc_mass_we(k,jy,j)+xmassperparticle/2.)/ &
                 xmassperparticle)
            acc_mass_we(k,jy,j)=acc_mass_we(k,jy,j)- &
                 real(mmass)*xmassperparticle
260
          else
261
            mmass=0
262
263
          endif

264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
          do m=1,mmass
            do ipart=minpart,maxpart

! If a vacant storage space is found, attribute everything to this array element
! TODO: for the MPI version this test can be removed, as all 
!       elements in _tmp arrays are initialized to zero
!*****************************************************************************

              if (itra1_tmp(ipart).ne.itime) then

! Assign particle positions
!**************************

                xtra1_tmp(ipart)=real(nx_we(k))
                if (jy.eq.ny_sn(1)) then
                  ytra1_tmp(ipart)=real(jy)+0.5*ran1(idummy)
                else if (jy.eq.ny_sn(2)) then
                  ytra1_tmp(ipart)=real(jy)-0.5*ran1(idummy)
                else
                  ytra1_tmp(ipart)=real(jy)+(ran1(idummy)-.5)
284
                endif
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
                if (j.eq.1) then
                  ztra1_tmp(ipart)=zcolumn_we(k,jy,1)+(zcolumn_we(k,jy,2)- &
                       zcolumn_we(k,jy,1))/4.
                else if (j.eq.numcolumn_we(k,jy)) then
                  ztra1_tmp(ipart)=(2.*zcolumn_we(k,jy,j)+ &
                       zcolumn_we(k,jy,j-1)+height(nz))/4.
                else
                  ztra1_tmp(ipart)=zcolumn_we(k,jy,j-1)+ran1(idummy)* &
                       (zcolumn_we(k,jy,j+1)-zcolumn_we(k,jy,j-1))
                endif

! Interpolate PV to the particle position
!****************************************
                ixm=int(xtra1_tmp(ipart))
                jym=int(ytra1_tmp(ipart))
                ixp=ixm+1
                jyp=jym+1
                ddx=xtra1_tmp(ipart)-real(ixm)
                ddy=ytra1_tmp(ipart)-real(jym)
                rddx=1.-ddx
                rddy=1.-ddy
                p1=rddx*rddy
                p2=ddx*rddy
                p3=rddx*ddy
                p4=ddx*ddy
                do i=2,nz
                  if (height(i).gt.ztra1_tmp(ipart)) then
                    indzm=i-1
                    indzp=i
                    goto 26
                  endif
316
                end do
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
26              continue
                dz1=ztra1_tmp(ipart)-height(indzm)
                dz2=height(indzp)-ztra1_tmp(ipart)
                dz=1./(dz1+dz2)
                do mm=1,2
                  indexh=memind(mm)
                  do in=1,2
                    indzh=indzm+in-1
                    y1(in)=p1*pv(ixm,jym,indzh,indexh) &
                         +p2*pv(ixp,jym,indzh,indexh) &
                         +p3*pv(ixm,jyp,indzh,indexh) &
                         +p4*pv(ixp,jyp,indzh,indexh)
                  end do
                  yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
                end do
                pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
                ylat=ylat0+ytra1_tmp(ipart)*dy
                if (ylat.lt.0.) pvpart=-1.*pvpart


! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
!*****************************************************************************

                if (((ztra1_tmp(ipart).gt.3000.).and. &
                     (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
                  nclass_tmp(ipart)=min(int(ran1(idummy)* &
                       real(nclassunc))+1,nclassunc)
                  numactiveparticles=numactiveparticles+1
                  numparticlecount=numparticlecount+1
                  npoint_tmp(ipart)=numparticlecount
                  idt_tmp(ipart)=mintime
                  itra1_tmp(ipart)=itime
                  itramem_tmp(ipart)=itra1_tmp(ipart)
                  itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit
                  xmass1_tmp(ipart,1)=xmassperparticle
                  if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= &
                       xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9
                else
                  goto 71
                endif
357
358


359
360
! Increase numpart, if necessary
!*******************************
361

362
363
364
365
366
367
368
369
                numpart_total=max(numpart_total,ipart)
                goto 73      ! Storage space has been found, stop searching
              endif
            end do
            if (ipart.gt.tmp_size) &
                 stop 'boundcond_domainfill_mpi.f90: too many particles required'
73          minpart=ipart+1
71          continue
370
371
372
          end do


373
        end do
374
375
376
377
      end do
    end do


378
379
380
!*****************************************
! Southern and northern boundary condition
!*****************************************
381

382
383
! Loop from west to east
!***********************
384

385
    do ix=nx_we(1),nx_we(2)
386

387
388
! Loop over southern (index 1) and northern (index 2) boundary
!*************************************************************
389

390
391
392
      do k=1,2
        ylat=ylat0+real(ny_sn(k))*dy
        cosfact=cos(ylat*pi180)
393

394
395
! Loop over all release locations in a column
!********************************************
396

397
        do j=1,numcolumn_sn(k,ix)
398

399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
! Determine, for each release location, the area of the corresponding boundary
!*****************************************************************************

          if (j.eq.1) then
            deltaz=(zcolumn_sn(k,ix,2)+zcolumn_sn(k,ix,1))/2.
          else if (j.eq.numcolumn_sn(k,ix)) then
!        deltaz=height(nz)-(zcolumn_sn(k,ix,j-1)+
!    +        zcolumn_sn(k,ix,j))/2.
! In order to avoid taking a very high column for very many particles,
! use the deltaz from one particle below instead
            deltaz=(zcolumn_sn(k,ix,j)-zcolumn_sn(k,ix,j-2))/2.
          else
            deltaz=(zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))/2.
          endif
          if ((ix.eq.nx_we(1)).or.(ix.eq.nx_we(2))) then
            boundarea=deltaz*111198.5/2.*cosfact*dx
          else
            boundarea=deltaz*111198.5*cosfact*dx
          endif
418
419


420
421
! Interpolate the wind velocity and density to the release location
!******************************************************************
422

423
424
! Determine the model level below the release position
!*****************************************************
425

426
427
428
429
430
431
432
433
          do i=2,nz
            if (height(i).gt.zcolumn_sn(k,ix,j)) then
              indz=i-1
              indzp=i
              goto 16
            endif
          end do
16        continue
434

435
436
! Vertical distance to the level below and above current position
!****************************************************************
437

438
439
440
          dz1=zcolumn_sn(k,ix,j)-height(indz)
          dz2=height(indzp)-zcolumn_sn(k,ix,j)
          dz=1./(dz1+dz2)
441

442
443
! Vertical and temporal interpolation
!************************************
444

445
446
447
448
449
450
451
          do m=1,2
            indexh=memind(m)
            do in=1,2
              indzh=indz+in-1
              windl(in)=vv(ix,ny_sn(k),indzh,indexh)
              rhol(in)=rho(ix,ny_sn(k),indzh,indexh)
            end do
452

453
454
            windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
            rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
455
456
          end do

457
458
          windx=(windhl(1)*dt2+windhl(2)*dt1)*dtt
          rhox=(rhohl(1)*dt2+rhohl(2)*dt1)*dtt
459

460
461
! Calculate mass flux
!********************
462

463
          fluxofmass=windx*rhox*boundarea*real(lsynctime)
464

465
466
467
! If the mass flux is directed into the domain, add it to previous mass fluxes;
! if it is out of the domain, set accumulated mass flux to zero
!******************************************************************************
468

469
470
471
472
473
474
          if (k.eq.1) then
            if (fluxofmass.ge.0.) then
              acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+fluxofmass
            else
              acc_mass_sn(k,ix,j)=0.
            endif
475
          else
476
477
478
479
480
            if (fluxofmass.le.0.) then
              acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)+abs(fluxofmass)
            else
              acc_mass_sn(k,ix,j)=0.
            endif
481
          endif
482
483
484
485
486
487
488
489
490
491
492
493
          accmasst=accmasst+acc_mass_sn(k,ix,j)

! If the accumulated mass exceeds half the mass that each particle shall carry,
! one (or more) particle(s) is (are) released and the accumulated mass is
! reduced by the mass of this (these) particle(s)
!******************************************************************************

          if (acc_mass_sn(k,ix,j).ge.xmassperparticle/2.) then
            mmass=int((acc_mass_sn(k,ix,j)+xmassperparticle/2.)/ &
                 xmassperparticle)
            acc_mass_sn(k,ix,j)=acc_mass_sn(k,ix,j)- &
                 real(mmass)*xmassperparticle
494
          else
495
            mmass=0
496
497
          endif

498
499
500
501
502
          do m=1,mmass
            do ipart=minpart,maxpart

! If a vacant storage space is found, attribute everything to this array element
!*****************************************************************************
503

504
505
506
507
508
509
510
511
512
513
514
515
              if (itra1_tmp(ipart).ne.itime) then

! Assign particle positions
!**************************

                ytra1_tmp(ipart)=real(ny_sn(k))
                if (ix.eq.nx_we(1)) then
                  xtra1_tmp(ipart)=real(ix)+0.5*ran1(idummy)
                else if (ix.eq.nx_we(2)) then
                  xtra1_tmp(ipart)=real(ix)-0.5*ran1(idummy)
                else
                  xtra1_tmp(ipart)=real(ix)+(ran1(idummy)-.5)
516
                endif
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
                if (j.eq.1) then
                  ztra1_tmp(ipart)=zcolumn_sn(k,ix,1)+(zcolumn_sn(k,ix,2)- &
                       zcolumn_sn(k,ix,1))/4.
                else if (j.eq.numcolumn_sn(k,ix)) then
                  ztra1_tmp(ipart)=(2.*zcolumn_sn(k,ix,j)+ &
                       zcolumn_sn(k,ix,j-1)+height(nz))/4.
                else
                  ztra1_tmp(ipart)=zcolumn_sn(k,ix,j-1)+ran1(idummy)* &
                       (zcolumn_sn(k,ix,j+1)-zcolumn_sn(k,ix,j-1))
                endif


! Interpolate PV to the particle position
!****************************************
                ixm=int(xtra1_tmp(ipart))
                jym=int(ytra1_tmp(ipart))
                ixp=ixm+1
                jyp=jym+1
                ddx=xtra1_tmp(ipart)-real(ixm)
                ddy=ytra1_tmp(ipart)-real(jym)
                rddx=1.-ddx
                rddy=1.-ddy
                p1=rddx*rddy
                p2=ddx*rddy
                p3=rddx*ddy
                p4=ddx*ddy
                do i=2,nz
                  if (height(i).gt.ztra1_tmp(ipart)) then
                    indzm=i-1
                    indzp=i
                    goto 126
                  endif
549
                end do
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
126             continue
                dz1=ztra1_tmp(ipart)-height(indzm)
                dz2=height(indzp)-ztra1_tmp(ipart)
                dz=1./(dz1+dz2)
                do mm=1,2
                  indexh=memind(mm)
                  do in=1,2
                    indzh=indzm+in-1
                    y1(in)=p1*pv(ixm,jym,indzh,indexh) &
                         +p2*pv(ixp,jym,indzh,indexh) &
                         +p3*pv(ixm,jyp,indzh,indexh) &
                         +p4*pv(ixp,jyp,indzh,indexh)
                  end do
                  yh1(mm)=(dz2*y1(1)+dz1*y1(2))*dz
                end do
                pvpart=(yh1(1)*dt2+yh1(2)*dt1)*dtt
                if (ylat.lt.0.) pvpart=-1.*pvpart


! For domain-filling option 2 (stratospheric O3), do the rest only in the stratosphere
!*****************************************************************************

                if (((ztra1_tmp(ipart).gt.3000.).and. &
                     (pvpart.gt.pvcrit)).or.(mdomainfill.eq.1)) then
                  nclass_tmp(ipart)=min(int(ran1(idummy)* &
                       real(nclassunc))+1,nclassunc)
                  numactiveparticles=numactiveparticles+1
                  numparticlecount=numparticlecount+1
                  npoint_tmp(ipart)=numparticlecount
                  idt_tmp(ipart)=mintime
                  itra1_tmp(ipart)=itime
                  itramem_tmp(ipart)=itra1_tmp(ipart)
                  itrasplit_tmp(ipart)=itra1_tmp(ipart)+ldirect*itsplit
                  xmass1_tmp(ipart,1)=xmassperparticle
                  if (mdomainfill.eq.2) xmass1_tmp(ipart,1)= &
                       xmass1_tmp(ipart,1)*pvpart*48./29.*ozonescale/10.**9
                else
                  goto 171
                endif
589
590


591
592
593
594
595
596
597
598
599
600
! Increase numpart, if necessary
!*******************************
                numpart_total=max(numpart_total,ipart)
                goto 173      ! Storage space has been found, stop searching
              endif
            end do
            if (ipart.gt.tmp_size) &
                 stop 'boundcond_domainfill.f: too many particles required'
173         minpart=ipart+1
171         continue
601
602
603
          end do


604
        end do
605
606
      end do
    end do
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625


! xm=0.
! do i=1,numpart_total
!   if (itra1_tmp(i).eq.itime) xm=xm+xmass1(i,1)
! end do

!write(*,*) itime,numactiveparticles,numparticlecount,numpart,
!    +xm,accmasst,xm+accmasst

  end if ! if lroot

! Distribute the number of particles to be released
! *************************************************
  call MPI_Bcast(numpart_total, 1, MPI_INTEGER, id_root, mp_comm_used, mp_ierr)

  do i=0, mp_partgroup_np-1
    numrel_mpi(i) = numpart_total/mp_partgroup_np
    if (i.lt.mod(numpart_total,mp_partgroup_np)) numrel_mpi(i) = numrel_mpi(i) + 1
626
627
  end do

628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
! Allocate temporary arrays for receiving processes
  if (.not.lroot) then 
    allocate(itra1_tmp(numrel_mpi(mp_partid)),&
         & npoint_tmp(numrel_mpi(mp_partid)),&
         & nclass_tmp(numrel_mpi(mp_partid)),&
         & idt_tmp(numrel_mpi(mp_partid)),&
         & itramem_tmp(numrel_mpi(mp_partid)),&
         & itrasplit_tmp(numrel_mpi(mp_partid)),&
         & xtra1_tmp(numrel_mpi(mp_partid)),&
         & ytra1_tmp(numrel_mpi(mp_partid)),&
         & ztra1_tmp(numrel_mpi(mp_partid)),&
         & xmass1_tmp(numrel_mpi(mp_partid),maxspec))
      
! Initialize all particles as non-existent    
    itra1_tmp(:)=-999999999
  end if
644

645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
! Distribute particles
! Keep track of released particles so far
  rel_counter = 0
  mtag = 1000

  do i=0, mp_partgroup_np-1

! For root process, nothing to do except update release count
    if (i.eq.0) then 
      rel_counter = rel_counter + numrel_mpi(i)
      cycle
    end if
      
! Send particles from root to non-root processes
    if (lroot.and.numrel_mpi(i).gt.0) then

      call MPI_SEND(nclass_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),MPI_INTEGER,i,mtag+1*i,mp_comm_used,mp_ierr)

      call MPI_SEND(npoint_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),MPI_INTEGER,i,mtag+2*i,mp_comm_used,mp_ierr)

      call MPI_SEND(itra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),MPI_INTEGER,i,mtag+3*i,mp_comm_used,mp_ierr)

      call MPI_SEND(idt_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),MPI_INTEGER,i,mtag+4*i,mp_comm_used,mp_ierr)

      call MPI_SEND(itramem_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),MPI_INTEGER,i,mtag+5*i,mp_comm_used,mp_ierr)

      call MPI_SEND(itrasplit_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),MPI_INTEGER,i,mtag+6*i,mp_comm_used,mp_ierr)

      call MPI_SEND(xtra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),mp_dp,i,mtag+7*i,mp_comm_used,mp_ierr)

      call MPI_SEND(ytra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),mp_dp,i,mtag+8*i,mp_comm_used,mp_ierr)

      call MPI_SEND(ztra1_tmp(rel_counter+1:rel_counter+numrel_mpi(i)),&
           &numrel_mpi(i),mp_sp,i,mtag+9*i,mp_comm_used,mp_ierr)

      do j=1,nspec
        call MPI_SEND(xmass1_tmp(rel_counter+1:rel_counter+numrel_mpi(i),j),&
             &numrel_mpi(i),mp_sp,i,mtag+(9+j)*i,mp_comm_used,mp_ierr)
      end do

! Non-root processes issue receive requests
    else if (i.eq.mp_partid.and.numrel_mpi(i).gt.0) then 
      call MPI_RECV(nclass_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &MPI_INTEGER,id_root,mtag+1*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(npoint_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &MPI_INTEGER,id_root,mtag+2*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(itra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &MPI_INTEGER,id_root,mtag+3*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(idt_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &MPI_INTEGER,id_root,mtag+4*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(itramem_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &MPI_INTEGER,id_root,mtag+5*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(itrasplit_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &MPI_INTEGER,id_root,mtag+6*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(xtra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &mp_dp,id_root,mtag+7*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(ytra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &mp_dp,id_root,mtag+8*i,mp_comm_used,mp_status,mp_ierr)

      call MPI_RECV(ztra1_tmp(1:numrel_mpi(i)),numrel_mpi(i),&
           &mp_sp,id_root,mtag+9*i,mp_comm_used,mp_status,mp_ierr)

      do j=1,nspec
        call MPI_RECV(xmass1_tmp(1:numrel_mpi(i),j),numrel_mpi(i),&
             &mp_sp,id_root,mtag+(9+j)*i,mp_comm_used,mp_status,mp_ierr)

      end do
    end if
    rel_counter = rel_counter + numrel_mpi(i)
729
730
  end do

731
732
733
! Find free storage space for the new particles.
! This section is independent of the redistribution scheme used
! ********************************************************************
734

735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
! Keep track of released particles so far
  minpart=1
    
! The algorithm should be correct also for root process
  do i=1, numrel_mpi(mp_partid)
    do ipart=minpart, maxpart
      if (itra1(ipart).ne.itime) then
        itra1(ipart) = itra1_tmp(i)
        npoint(ipart) = npoint_tmp(i)
        nclass(ipart) = nclass_tmp(i)
        idt(ipart) = idt_tmp(i)
        itramem(ipart) = itramem_tmp(i)
        itrasplit(ipart) = itrasplit_tmp(i)
        xtra1(ipart) = xtra1_tmp(i)
        ytra1(ipart) = ytra1_tmp(i)
        ztra1(ipart) = ztra1_tmp(i)
        xmass1(ipart,:) = xmass1_tmp(i,:)
! Increase numpart, if necessary
        numpart=max(numpart,ipart)
        goto 200 ! Storage space has been found, stop searching
      end if
    end do
200 minpart=ipart+1
  end do
759

760
761
762
! If particles shall be dumped, then accumulated masses at the domain boundaries
! must be dumped, too, to be used for later runs
!*****************************************************************************
763
764

  if ((ipout.gt.0).and.(itime.eq.loutend)) then
765
    if (lroot) then
766
      call mpif_mtime('iotime',0)
767
768
769
770
771
      open(unitboundcond,file=path(2)(1:length(2))//'boundcond.bin', &
           form='unformatted')
      write(unitboundcond) numcolumn_we,numcolumn_sn, &
           zcolumn_we,zcolumn_sn,acc_mass_we,acc_mass_sn
      close(unitboundcond)
772
      call mpif_mtime('iotime',1)
773
    end if
774
775
  endif

776
777
778
779
780
781
! Deallocate temporary arrays 
  deallocate(itra1_tmp,npoint_tmp,nclass_tmp,idt_tmp,itramem_tmp,itrasplit_tmp,&
       & xtra1_tmp,ytra1_tmp,ztra1_tmp,xmass1_tmp,numrel_mpi)
! numactiveparticles_mpi


782
end subroutine boundcond_domainfill