boundcond_domainfill_mpi.f90 27.9 KB
Newer Older
1
subroutine boundcond_domainfill(itime,loutend)
2
3
4
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
!                                  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)     
!*****************************************************************************
35
36
37
38
39
40
41
42
43
44
45
46

  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
47
48
  integer :: numactiveparticles, numpart_total, rel_counter
  integer,allocatable,dimension(:) ::  numrel_mpi !, numactiveparticles_mpi
49
50
51
52
53
54
55
56
57
58

  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
59
  integer :: mtag
60
  logical :: first_call=.true.
61
62
63
64
! Sizes of temporary arrays are maxpartfract*maxpart. Increase maxpartfract if 
! needed.  
  real,parameter :: maxpartfract=0.1
  integer :: tmp_size = int(maxpartfract*maxpart)
65

66
! Use different seed for each process
67
68
69
70
71
72
  if (first_call) then
    idummy=idummy+mp_seed
    first_call=.false.
  end if


73
74
! If domain-filling is global, no boundary conditions are needed
!***************************************************************
75
76
77
78
79

  if (gdomainfill) return

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

84
85
86
87
88
89
! 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
!********************************************************************
90
91
92
93
94
95
96
97
98
99
100
101

  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
102
!  numactiveparticles_mpi(mp_partid) = numactiveparticles
103
104


105
106
107
! 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)
108
109


110
111
! Total number of new releases
  numpart_total = 0
112
113


114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
! 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
!*****************************************************
132

133
134
135
    dt1=real(itime-memtime(1))
    dt2=real(memtime(2)-itime)
    dtt=1./(dt1+dt2)
136

137
138
! Initialize auxiliary variable used to search for vacant storage space
!**********************************************************************
139

140
    minpart=1
141

142
143
144
!***************************************
! Western and eastern boundary condition
!***************************************
145

146
147
! Loop from south to north
!*************************
148

149
    do jy=ny_sn(1),ny_sn(2)
150

151
152
! Loop over western (index 1) and eastern (index 2) boundary
!***********************************************************
153

154
      do k=1,2
155

156
157
! Loop over all release locations in a column
!********************************************
158

159
        do j=1,numcolumn_we(k,jy)
160

161
162
! Determine, for each release location, the area of the corresponding boundary
!*****************************************************************************
163

164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
          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
179
180
181
          endif


182
183
! Interpolate the wind velocity and density to the release location
!******************************************************************
184

185
186
! Determine the model level below the release position
!*****************************************************
187

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

197
198
! Vertical distance to the level below and above current position
!****************************************************************
199

200
201
202
          dz1=zcolumn_we(k,jy,j)-height(indz)
          dz2=height(indzp)-zcolumn_we(k,jy,j)
          dz=1./(dz1+dz2)
203

204
205
! Vertical and temporal interpolation
!************************************
206

207
208
209
210
211
212
213
          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
214

215
216
217
218
219
220
221
222
223
224
225
            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)
226
227


228
229
230
231
232
233
234
235
236
237
! 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
238
          else
239
240
241
242
243
            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
244
          endif
245
246
247
248
249
250
251
252
253
254
255
256
          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
257
          else
258
            mmass=0
259
260
          endif

261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
          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)
281
                endif
282
283
284
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
                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
313
                end do
314
315
316
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
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
354
355


356
357
! Increase numpart, if necessary
!*******************************
358

359
360
361
362
363
364
365
366
                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
367
368
369
          end do


370
        end do
371
372
373
374
      end do
    end do


375
376
377
!*****************************************
! Southern and northern boundary condition
!*****************************************
378

379
380
! Loop from west to east
!***********************
381

382
    do ix=nx_we(1),nx_we(2)
383

384
385
! Loop over southern (index 1) and northern (index 2) boundary
!*************************************************************
386

387
388
389
      do k=1,2
        ylat=ylat0+real(ny_sn(k))*dy
        cosfact=cos(ylat*pi180)
390

391
392
! Loop over all release locations in a column
!********************************************
393

394
        do j=1,numcolumn_sn(k,ix)
395

396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
! 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
415
416


417
418
! Interpolate the wind velocity and density to the release location
!******************************************************************
419

420
421
! Determine the model level below the release position
!*****************************************************
422

423
424
425
426
427
428
429
430
          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
431

432
433
! Vertical distance to the level below and above current position
!****************************************************************
434

435
436
437
          dz1=zcolumn_sn(k,ix,j)-height(indz)
          dz2=height(indzp)-zcolumn_sn(k,ix,j)
          dz=1./(dz1+dz2)
438

439
440
! Vertical and temporal interpolation
!************************************
441

442
443
444
445
446
447
448
          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
449

450
451
            windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
            rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
452
453
          end do

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

457
458
! Calculate mass flux
!********************
459

460
          fluxofmass=windx*rhox*boundarea*real(lsynctime)
461

462
463
464
! 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
!******************************************************************************
465

466
467
468
469
470
471
          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
472
          else
473
474
475
476
477
            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
478
          endif
479
480
481
482
483
484
485
486
487
488
489
490
          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
491
          else
492
            mmass=0
493
494
          endif

495
496
497
498
499
          do m=1,mmass
            do ipart=minpart,maxpart

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

501
502
503
504
505
506
507
508
509
510
511
512
              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)
513
                endif
514
515
516
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
                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
546
                end do
547
548
549
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
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
586
587


588
589
590
591
592
593
594
595
596
597
! 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
598
599
600
          end do


601
        end do
602
603
      end do
    end do
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622


! 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
623
624
  end do

625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
! 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
641

642
643
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
! 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)
726
727
  end do

728
729
730
! Find free storage space for the new particles.
! This section is independent of the redistribution scheme used
! ********************************************************************
731

732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
! 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
756

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

  if ((ipout.gt.0).and.(itime.eq.loutend)) then
762
    if (lroot) then
763
      call mpif_mtime('iotime',0)
764
765
766
767
768
      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)
769
      call mpif_mtime('iotime',1)
770
    end if
771
772
  endif

773
774
775
776
777
778
! 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


779
end subroutine boundcond_domainfill