boundcond_domainfill_mpi.f90 29.3 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

subroutine boundcond_domainfill(itime,loutend)
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
!                                  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)     
!*****************************************************************************
56
57
58
59
60
61
62
63
64
65
66
67

  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
68
69
  integer :: numactiveparticles, numpart_total, rel_counter
  integer,allocatable,dimension(:) ::  numrel_mpi !, numactiveparticles_mpi
70
71
72
73
74
75
76
77
78
79

  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
80
  integer :: mtag
81
  logical :: first_call=.true.
82
83
84
85
! Sizes of temporary arrays are maxpartfract*maxpart. Increase maxpartfract if 
! needed.  
  real,parameter :: maxpartfract=0.1
  integer :: tmp_size = int(maxpartfract*maxpart)
86

87
! Use different seed for each process
88
89
90
91
92
93
  if (first_call) then
    idummy=idummy+mp_seed
    first_call=.false.
  end if


94
95
! If domain-filling is global, no boundary conditions are needed
!***************************************************************
96
97
98
99
100

  if (gdomainfill) return

  accmasst=0.
  numactiveparticles=0
101
102
103
! Keep track of active particles on each process
  allocate(numrel_mpi(0:mp_partgroup_np-1))
! numactiveparticles_mpi(0:mp_partgroup_np-1)
104

105
106
107
108
109
110
! 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
!********************************************************************
111
112
113
114
115
116
117
118
119
120
121
122

  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
123
!  numactiveparticles_mpi(mp_partid) = numactiveparticles
124
125


126
127
128
! 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)
129
130


131
132
! Total number of new releases
  numpart_total = 0
133
134


135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
! 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
!*****************************************************
153

154
155
156
    dt1=real(itime-memtime(1))
    dt2=real(memtime(2)-itime)
    dtt=1./(dt1+dt2)
157

158
159
! Initialize auxiliary variable used to search for vacant storage space
!**********************************************************************
160

161
    minpart=1
162

163
164
165
!***************************************
! Western and eastern boundary condition
!***************************************
166

167
168
! Loop from south to north
!*************************
169

170
    do jy=ny_sn(1),ny_sn(2)
171

172
173
! Loop over western (index 1) and eastern (index 2) boundary
!***********************************************************
174

175
      do k=1,2
176

177
178
! Loop over all release locations in a column
!********************************************
179

180
        do j=1,numcolumn_we(k,jy)
181

182
183
! Determine, for each release location, the area of the corresponding boundary
!*****************************************************************************
184

185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
          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
200
201
202
          endif


203
204
! Interpolate the wind velocity and density to the release location
!******************************************************************
205

206
207
! Determine the model level below the release position
!*****************************************************
208

209
210
211
212
213
214
          do i=2,nz
            if (height(i).gt.zcolumn_we(k,jy,j)) then
              indz=i-1
              indzp=i
              goto 6
            endif
215
          end do
216
6         continue
217

218
219
! Vertical distance to the level below and above current position
!****************************************************************
220

221
222
223
          dz1=zcolumn_we(k,jy,j)-height(indz)
          dz2=height(indzp)-zcolumn_we(k,jy,j)
          dz=1./(dz1+dz2)
224

225
226
! Vertical and temporal interpolation
!************************************
227

228
229
230
231
232
233
234
          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
235

236
237
238
239
240
241
242
243
244
245
246
            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)
247
248


249
250
251
252
253
254
255
256
257
258
! 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
259
          else
260
261
262
263
264
            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
265
          endif
266
267
268
269
270
271
272
273
274
275
276
277
          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
278
          else
279
            mmass=0
280
281
          endif

282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
          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)
302
                endif
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
                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
334
                end do
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
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
375
376


377
378
! Increase numpart, if necessary
!*******************************
379

380
381
382
383
384
385
386
387
                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
388
389
390
          end do


391
        end do
392
393
394
395
      end do
    end do


396
397
398
!*****************************************
! Southern and northern boundary condition
!*****************************************
399

400
401
! Loop from west to east
!***********************
402

403
    do ix=nx_we(1),nx_we(2)
404

405
406
! Loop over southern (index 1) and northern (index 2) boundary
!*************************************************************
407

408
409
410
      do k=1,2
        ylat=ylat0+real(ny_sn(k))*dy
        cosfact=cos(ylat*pi180)
411

412
413
! Loop over all release locations in a column
!********************************************
414

415
        do j=1,numcolumn_sn(k,ix)
416

417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
! 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
436
437


438
439
! Interpolate the wind velocity and density to the release location
!******************************************************************
440

441
442
! Determine the model level below the release position
!*****************************************************
443

444
445
446
447
448
449
450
451
          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
452

453
454
! Vertical distance to the level below and above current position
!****************************************************************
455

456
457
458
          dz1=zcolumn_sn(k,ix,j)-height(indz)
          dz2=height(indzp)-zcolumn_sn(k,ix,j)
          dz=1./(dz1+dz2)
459

460
461
! Vertical and temporal interpolation
!************************************
462

463
464
465
466
467
468
469
          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
470

471
472
            windhl(m)=(dz2*windl(1)+dz1*windl(2))*dz
            rhohl(m)=(dz2*rhol(1)+dz1*rhol(2))*dz
473
474
          end do

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

478
479
! Calculate mass flux
!********************
480

481
          fluxofmass=windx*rhox*boundarea*real(lsynctime)
482

483
484
485
! 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
!******************************************************************************
486

487
488
489
490
491
492
          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
493
          else
494
495
496
497
498
            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
499
          endif
500
501
502
503
504
505
506
507
508
509
510
511
          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
512
          else
513
            mmass=0
514
515
          endif

516
517
518
519
520
          do m=1,mmass
            do ipart=minpart,maxpart

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

522
523
524
525
526
527
528
529
530
531
532
533
              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)
534
                endif
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
                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
567
                end do
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
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
607
608


609
610
611
612
613
614
615
616
617
618
! 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
619
620
621
          end do


622
        end do
623
624
      end do
    end do
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643


! 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
644
645
  end do

646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
! 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
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
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
! 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)
747
748
  end do

749
750
751
! Find free storage space for the new particles.
! This section is independent of the redistribution scheme used
! ********************************************************************
752

753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
! 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
777

778
779
780
! If particles shall be dumped, then accumulated masses at the domain boundaries
! must be dumped, too, to be used for later runs
!*****************************************************************************
781
782

  if ((ipout.gt.0).and.(itime.eq.loutend)) then
783
    if (lroot) then
784
      call mpif_mtime('iotime',0)
785
786
787
788
789
      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)
790
      call mpif_mtime('iotime',1)
791
    end if
792
793
  endif

794
795
796
797
798
799
! 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


800
end subroutine boundcond_domainfill