source: palm/trunk/SOURCE/data_output_particle_mod.f90 @ 4782

Last change on this file since 4782 was 4779, checked in by suehring, 4 years ago

Avoid overlong lines and unused variables

  • Property svn:keywords set to Id
File size: 82.9 KB
Line 
1!> @file data_output_particle_mod.f90
2!--------------------------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
6! Public License as published by the Free Software Foundation, either version 3 of the License, or
7! (at your option) any later version.
8!
9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
11! Public License for more details.
12!
13! You should have received a copy of the GNU General Public License along with PALM. If not, see
14! <http://www.gnu.org/licenses/>.
15!
16! Copyright 1997-2020 Leibniz Universitaet Hannover
17!--------------------------------------------------------------------------------------------------!
18!
19! Current revisions:
20! -----------------
21!
22!
23! Former revisions:
24! -----------------
25! $Id: data_output_particle_mod.f90 4779 2020-11-09 17:45:22Z gronemeier $
26! Avoid overlong lines and unused variables
27!
28! 4778 2020-11-09 13:40:05Z raasch
29! Initial implementation (K. Ketelsen)
30!
31!
32! Description:
33! ------------
34!> Output of particle time series
35!--------------------------------------------------------------------------------------------------!
36 MODULE data_output_particle_mod
37
38#if defined( __parallel )
39   USE MPI
40#endif
41
42#if defined( __netcdf4 )
43   USE NETCDF
44#endif
45
46   USE, INTRINSIC ::  ISO_C_BINDING
47
48   USE kinds,                                                                                     &
49      ONLY: wp, iwp, sp, dp, idp, isp
50
51   USE indices,                                                                                   &
52       ONLY:  nbgp, nnx, nny, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz, nzb, nzt
53
54   USE control_parameters,                                                                        &
55       ONLY:   end_time, simulated_time, dt_dopts
56
57   USE pegrid,                                                                                    &
58       ONLY:  comm1dx, comm1dy, comm2d, myid, myidx, myidy, npex, npey, numprocs, pdims
59
60   USE particle_attributes,                                                                       &
61       ONLY: particle_type, particles, grid_particle_def, grid_particles, prt_count,              &
62             unlimited_dimension, pts_id_file, data_output_pts, pts_increment, pts_percentage, &
63             oversize, number_of_output_particles
64
65   USE shared_memory_io_mod,                                                                      &
66       ONLY: sm_class
67
68   USE data_output_netcdf4_module,                                                                   &
69       ONLY: netcdf4_init_dimension, netcdf4_get_error_message, netcdf4_stop_file_header_definition, &
70             netcdf4_init_module, netcdf4_init_variable, netcdf4_finalize,                           &
71             netcdf4_open_file, netcdf4_write_attribute, netcdf4_write_variable
72
73    USE cpulog,                                                                &
74        ONLY:  cpu_log, log_point_s
75
76   IMPLICIT NONE
77
78   PRIVATE
79   SAVE
80
81
82   LOGICAL, PUBLIC      :: dop_active = .FALSE.
83
84
85!  Variables for restart
86
87   INTEGER(iwp), PUBLIC    :: dop_prt_axis_dimension
88   INTEGER(iwp), PUBLIC    :: dop_last_active_particle
89
90
91
92
93!kk Private module variables can be declared inside the submodule
94
95   INTEGER,PARAMETER    :: MAX_NR_VARIABLES    = 128
96   INTEGER,PARAMETER    :: NR_FIXED_VARIABLES  = 16
97
98   CHARACTER(LEN=32)                         :: file_name            !< Name of NetCDF file
99   INTEGER(iwp)                              :: nr_time_values       !< Number of values on time axis
100   REAL(sp),ALLOCATABLE,DIMENSION(:)         :: time_axis_values     !< time axis Values
101   INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)   :: rma_particles        !< Start address and number of remote particles
102   INTEGER(iwp)                              :: nr_particles_PE      !< Number of particles assigned for output on this thread
103   INTEGER(iwp)                              :: nr_particles_out     !< total number od particles assigned for output
104
105   INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)   :: sh_indices                   !< Indices in shared memory group
106   INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)   :: io_indices                   !< Indices on IO processes
107   INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)   :: mo_indices                   !< Indices for model communicator
108   INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)   :: remote_nr_particles
109   INTEGER(iwp)                              :: start_local_numbering        !< start increment 1 numbering on this thread
110   INTEGER(iwp)                              :: end_local_numbering
111   INTEGER(iwp)                              :: pe_start_index               !< start index of the output area on this thread
112   INTEGER(iwp)                              :: pe_end_index
113   INTEGER(iwp)                              :: io_start_index               !< start index of the output area on IO thread
114   INTEGER(iwp)                              :: io_end_index
115   INTEGER(iwp)                              :: nr_particles_rest            !< Numbere of rest Particle (ireg. distribution)
116   LOGICAL                                   :: irregular_distribubtion      !< irregular distribution of output particlesexit
117
118   TYPE var_def
119      INTEGER(iwp)         :: var_id
120      CHARACTER(len=32)    :: name
121      CHARACTER(len=32)    :: units
122      LOGICAL              :: is_integer = .FALSE.
123   END TYPE var_def
124
125   TYPE(var_def),DIMENSION(MAX_NR_VARIABLES)    :: variables
126   TYPE(var_def),DIMENSION(NR_FIXED_VARIABLES)  :: fix_variables
127
128   TYPE(sm_class)  :: part_io                 !< manage communicator for particle IO
129
130!
131!  NetCDF
132
133   INTEGER(iwp)                              :: file_id = -1         !< id of Netcdf file
134   INTEGER(iwp)                              :: nr_fix_variables     !< Number of fixed variables scheduled for output
135   INTEGER(iwp)                              :: nr_variables         !< Number of variables  scheduled for output
136
137   TYPE dimension_id
138      INTEGER(iwp)         :: prt
139      INTEGER(iwp)         :: time
140   END TYPE dimension_id
141
142   TYPE variable_id
143      INTEGER(iwp)         :: prt
144      INTEGER(iwp)         :: time
145   END TYPE variable_id
146
147   TYPE(dimension_id)                        :: did
148   TYPE(variable_id)                         :: var_id
149
150!  shared memory buffer
151   INTEGER(iwp)    :: win_prt_i = -1                                 ! integer MPI shared Memory window
152   INTEGER(iwp)    :: win_prt_r = -1                                 !< real  MPI shared Memory window
153   INTEGER(iwp), POINTER, CONTIGUOUS, DIMENSION(:) :: out_buf_i      !< integer output buffer
154   REAL(sp), POINTER, CONTIGUOUS, DIMENSION(:)     :: out_buf_r      !< real output buffer
155
156!
157!  Particle list in file
158
159   LOGICAL         :: part_list_in_file
160   INTEGER(idp),ALLOCATABLE, DIMENSION(:)   :: part_id_list_file
161!
162!  RMA window
163#if defined( __parallel )
164   INTEGER(iwp)    ::  win_rma_buf_i
165   INTEGER(iwp)    ::  win_rma_buf_r
166#endif
167
168   INTEGER(iwp),POINTER,DIMENSION(:)         :: transfer_buffer_i    !< rma window to provide data, which can be fetch via MPI_Get
169   REAL(sp),POINTER,DIMENSION(:)             :: transfer_buffer_r    !< Same for REAL
170   INTEGER(iwp),ALLOCATABLE,DIMENSION(:)     :: remote_indices       !< particle nubmer of the remodte partices, used as indices in the output array
171   INTEGER(iwp)                              :: initial_number_of_active_particles
172
173!  Public subroutine Interface
174
175   INTERFACE dop_init
176      MODULE PROCEDURE dop_init
177   END INTERFACE dop_init
178
179   INTERFACE dop_output_tseries
180      MODULE PROCEDURE dop_output_tseries
181   END INTERFACE dop_output_tseries
182
183   INTERFACE dop_finalize
184      MODULE PROCEDURE dop_finalize
185   END INTERFACE  dop_finalize
186
187#if defined( __parallel )
188   INTERFACE dop_alloc_rma_mem
189      MODULE PROCEDURE dop_alloc_rma_mem_i1
190      MODULE PROCEDURE dop_alloc_rma_mem_r1
191   END INTERFACE dop_alloc_rma_mem
192#endif
193
194   PUBLIC dop_init, dop_output_tseries, dop_finalize
195#if defined( __parallel )
196   PUBLIC dop_alloc_rma_mem       ! Must be PUBLIC on NEC, although it is only used in Submodule
197#endif
198
199
200
201 CONTAINS
202
203   SUBROUTINE dop_init (read_restart)
204      IMPLICIT NONE
205      LOGICAL,INTENT(IN)    :: read_restart
206
207      INTEGER(iwp) :: i                           !<
208      INTEGER(iwp) :: nr_particles_local          !< total number of particles scheduled for output on this thread
209#if defined( __parallel )
210      INTEGER(iwp) :: ierr                        !< MPI error code
211#endif
212      INTEGER(idp) :: nr_particles_8              !< Total number of particles in 64 bit
213      REAL(dp)     :: xnr_part                    !< Must be 64 Bit REAL
214      INTEGER(iwp) :: nr_local_last_pe               !< Number of output particles on myid == numprocs-2
215
216      INTEGER(iwp),DIMENSION(0:numprocs-1)     :: nr_particles_all_s
217      INTEGER(iwp),DIMENSION(0:numprocs-1)     :: nr_particles_all_r
218      INTEGER(idp),DIMENSION(0:numprocs-1)     :: nr_particles_all_8
219
220      INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)  :: sh_indices_s
221      INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)  :: io_indices_s
222      INTEGER(iwp),ALLOCATABLE,DIMENSION(:,:)  :: mo_indices_s
223
224      part_list_in_file = (pts_id_file /= ' ')
225
226      IF(.NOT. part_list_in_file)   THEN
227         IF(.NOT. read_restart)   THEN
228
229            CALL set_indef_particle_nr
230
231            CALL count_output_particles (nr_particles_local)
232
233            nr_particles_all_s = 0
234            nr_particles_all_s(myid) = nr_particles_local
235
236#if defined( __parallel )
237            CALL MPI_ALLREDUCE( nr_particles_all_s, nr_particles_all_r, SIZE( nr_particles_all_s ), MPI_INTEGER,    &
238               MPI_SUM, comm2d, ierr )
239#else
240            nr_particles_all_r = nr_particles_all_s
241#endif
242
243            initial_number_of_active_particles = SUM(nr_particles_all_r)
244
245            start_local_numbering = 1
246            end_local_numbering   = nr_particles_all_r(0)
247
248            IF ( myid > 0)   THEN
249               DO i=1,numprocs-1
250                  start_local_numbering = start_local_numbering+nr_particles_all_r(i-1)
251                  end_local_numbering   = start_local_numbering+nr_particles_all_r(i)-1
252                  IF(myid == i)  EXIT
253               ENDDO
254            END IF
255
256            nr_particles_all_8 = nr_particles_all_r              ! Use 64 INTEGER, maybe more than 2 GB particles
257            nr_particles_8 = SUM(nr_particles_all_8)
258            IF (nr_particles_8 > 2000000000_8)   THEN            ! This module works only with 32 Bit INTEGER
259               write(9,*) 'Number of particles too large ',nr_particles_8; flush(9)
260#if defined( __parallel )
261               CALL MPI_ABORT (MPI_COMM_WORLD, 1, ierr)
262#endif
263            ENDIF
264            nr_particles_out = nr_particles_8                    ! Total number of particles scheduled for output
265
266!
267!--         reserve space for additional particle
268            IF(number_of_output_particles > 0)   THEN
269               nr_particles_out = MAX(number_of_output_particles, nr_particles_out)
270            ELSE IF(oversize > 100.)   THEN
271               xnr_part = nr_particles_out*oversize/100.
272               nr_particles_out = xnr_part
273            ENDIF
274         ELSE
275            nr_particles_out = dop_prt_axis_dimension            ! Get Number from restart file
276            initial_number_of_active_particles = dop_last_active_particle
277         ENDIF
278
279      ELSE
280         CALL set_indef_particle_nr
281         CALL dop_read_output_particle_list (nr_particles_out)
282      ENDIF
283      dop_prt_axis_dimension = nr_particles_out
284!
285!--   The number of particles must be at least the number of MPI processes
286
287      nr_particles_out = MAX(nr_particles_out, numprocs)
288
289
290      nr_particles_PE = (nr_particles_out+numprocs-1)/numprocs   ! Number of paricles scheduled for ouput on this thread
291
292      pe_start_index = myid*nr_particles_PE+1                              !Output numberimng on this thread
293      pe_end_index   = MIN((myid+1)*nr_particles_PE,nr_particles_out)
294
295      irregular_distribubtion = .FALSE.
296#if defined( __parallel )
297!
298!--   In case of few particles, it can happen that not only the last thread gets fewer output particles.
299!--   In this case, the local number of particles on thread numprocs-1 will be < 1
300!--   If this happens, irregular distribution of output particles will be used
301
302      IF(myid == numprocs-1)   Then
303         nr_local_last_pe = pe_end_index-pe_start_index+1
304      ELSE
305         nr_local_last_pe = 0
306      ENDIF
307      CALL MPI_BCAST (nr_local_last_pe, 1, MPI_INTEGER, numprocs-1, comm2d, ierr)
308#else
309      nr_local_last_pe = nr_particles_PE
310#endif
311      IF(nr_local_last_pe < 1)   THEN
312         irregular_distribubtion = .TRUE.
313         CALL dop_setup_ireg_distribution
314      ENDIF
315
316
317      IF(.NOT. read_restart .AND. .NOT. part_list_in_file)   THEN
318         CALL set_particle_number
319      ENDIF
320
321      IF(part_list_in_file)   THEN
322         CALL dop_find_particle_in_outlist
323      ENDIF
324
325      CALL part_io%sm_init_data_output_particles ()
326
327      ALLOCATE (sh_indices_s(2,0:part_io%sh_npes-1))
328      ALLOCATE (sh_indices(2,0:part_io%sh_npes-1))
329
330
331#if defined( __parallel )
332      sh_indices_s = 0
333      sh_indices_s(1,part_io%sh_rank) = pe_start_index
334      sh_indices_s(2,part_io%sh_rank) = pe_end_index
335
336      CALL MPI_ALLREDUCE( sh_indices_s, sh_indices, 2*part_io%sh_npes, MPI_INTEGER,    &
337                                              MPI_SUM, part_io%comm_shared, ierr )
338
339      io_start_index = sh_indices(1,0)                         ! output numbering on actual IO thread
340      io_end_index   = sh_indices(2,part_io%sh_npes-1)
341#else
342      io_start_index = pe_start_index                          ! output numbering
343      io_end_index   = pe_end_index
344#endif
345
346
347#if defined( __parallel )
348      CALL MPI_BCAST (part_io%io_npes, 1, MPI_INTEGER, 0,  part_io%comm_shared, ierr)
349#endif
350      ALLOCATE (io_indices(2,0:part_io%io_npes-1))
351      IF (part_io%iam_io_pe)   THEN
352         ALLOCATE (io_indices_s(2,0:part_io%io_npes-1))
353
354
355         io_indices_s = 0
356         io_indices_s(1,part_io%io_rank) = io_start_index
357         io_indices_s(2,part_io%io_rank) = io_end_index
358
359#if defined( __parallel )
360         CALL MPI_ALLREDUCE( io_indices_s, io_indices, 2*part_io%io_npes, MPI_INTEGER,    &
361            MPI_SUM, part_io%comm_io, ierr )
362#else
363         io_indices = io_indices_s
364#endif
365
366      ENDIF
367
368#if defined( __parallel )
369      CALL MPI_BCAST (io_indices, size(io_indices), MPI_INTEGER, 0,  part_io%comm_shared, ierr)
370#endif
371
372      ALLOCATE (remote_nr_particles(2,0:numprocs-1))
373      ALLOCATE (rma_particles(2,0:numprocs-1))
374
375      ALLOCATE (mo_indices(2,0:numprocs-1))
376      ALLOCATE (mo_indices_s(2,0:numprocs-1))
377
378
379      mo_indices_s = 0
380      mo_indices_s(1,myid) = pe_start_index
381      mo_indices_s(2,myid) = pe_end_index
382
383#if defined( __parallel )
384      CALL MPI_ALLREDUCE( mo_indices_s, mo_indices, 2*numprocs, MPI_INTEGER, MPI_SUM, comm2d, ierr )
385#else
386      mo_indices = mo_indices_s
387#endif
388
389!--   Allocate output buffer
390
391#if defined( __parallel )
392      CALL part_io%sm_allocate_shared (out_buf_r, io_start_index, io_end_index, win_prt_r)
393      CALL part_io%sm_allocate_shared (out_buf_i, io_start_index, io_end_index, win_prt_i)
394#else
395      ALLOCATE(out_buf_r(io_start_index:io_end_index))
396      ALLOCATE(out_buf_i(io_start_index:io_end_index))
397#endif
398
399!--   NetCDF
400
401      CALL dop_netcdf_setup ()
402
403#if defined( __parallel )
404      CALL MPI_BCAST (nr_fix_variables, 1, MPI_INTEGER, 0,  part_io%comm_shared, ierr)
405      CALL MPI_BCAST (nr_variables, 1, MPI_INTEGER, 0,  part_io%comm_shared, ierr)
406#endif
407
408      CALL dop_count_remote_particles
409
410      CALL dop_write_fixed_variables
411
412      CALL deallocate_and_free
413
414      CALL dop_output_tseries
415
416      RETURN
417   CONTAINS
418      SUBROUTINE dop_setup_ireg_distribution
419         IMPLICIT NONE
420
421
422         nr_particles_PE   = (nr_particles_out)/numprocs   ! Number of paricles scheduled for ouput on this thread
423
424         nr_particles_rest = nr_particles_out - numprocs * nr_particles_PE
425
426         nr_particles_PE = nr_particles_PE+1
427         IF(myid < nr_particles_rest)   THEN
428            pe_start_index = myid*nr_particles_PE+1                              !Output numberimng on this thread
429            pe_end_index   = MIN((myid+1)*nr_particles_PE,nr_particles_out)
430         ELSE
431            pe_start_index = nr_particles_rest*(nr_particles_PE)+(myid-nr_particles_rest)*(nr_particles_PE-1)+1                              !Output numberimng on this thread
432            pe_end_index   = MIN(pe_start_index+nr_particles_PE-2,nr_particles_out)
433         ENDIF
434
435
436      END SUBROUTINE dop_setup_ireg_distribution
437
438   END SUBROUTINE dop_init
439!
440!  particle output on selected time steps
441
442   SUBROUTINE dop_output_tseries
443      IMPLICIT NONE
444
445      INTEGER(iwp)             :: i                           !<
446      INTEGER(iwp)             :: return_value                !< Return value data_output_netcdf4 .. routines
447#if defined( __parallel )
448      INTEGER(iwp)             :: ierr                        !< MPI error code
449#endif
450      INTEGER(iwp),SAVE        :: icount=0                    !< count output steps
451
452      INTEGER(iwp),DIMENSION(2)      :: bounds_origin, bounds_start, value_counts
453      REAl(wp),POINTER, CONTIGUOUS, DIMENSION(:)  :: my_time
454
455      icount = icount+1
456
457      CALL dop_delete_particle_number
458
459      IF(part_list_in_file)   THEN
460
461         CALL dop_find_particle_in_outlist
462      ELSE
463         CALL dop_newly_generated_particles
464      ENDIF
465
466      CALL dop_count_remote_particles
467
468      bounds_origin    = 1
469
470      bounds_start(1)  = io_start_index
471      bounds_start(2)  = icount
472
473      value_counts(1)  = io_end_index-io_start_index+1
474      value_counts(2)  = 1
475
476      DO i=1,nr_variables
477#if defined( __netcdf4 )
478         IF(variables(i)%is_integer)   THEN
479            out_buf_i = NF90_FILL_INT
480         ELSE
481            out_buf_r = NF90_FILL_REAL
482         ENDIF
483#endif
484
485         CALL cpu_log( log_point_s(99), 'dop_fill_out_buf', 'start' )
486         CALL dop_fill_out_buf (variables(i))
487         CALL cpu_log( log_point_s(99), 'dop_fill_out_buf', 'stop' )
488
489         CALL cpu_log( log_point_s(88), 'dop_get_remote_particle', 'start' )
490#if defined( __parallel )
491         CALL  dop_get_remote_particle (variables(i)%is_integer)
492#endif
493         CALL cpu_log( log_point_s(88), 'dop_get_remote_particle', 'stop' )
494
495         CALL cpu_log( log_point_s(89), 'particle NetCDF output', 'start' )
496         CALL  part_io%sm_node_barrier()
497         IF (part_io%iam_io_pe)   THEN
498            IF(variables(i)%is_integer)   THEN
499               CALL  netcdf4_write_variable('parallel', file_id, variables(i)%var_id, bounds_start,&
500                                             value_counts, bounds_origin,                          &
501                                             .FALSE., values_int32_1d=out_buf_i, return_value=return_value)
502            ELSE
503               CALL  netcdf4_write_variable('parallel', file_id, variables(i)%var_id, bounds_start,&
504                                             value_counts, bounds_origin,                          &
505                                            .FALSE., values_real32_1d=out_buf_r, return_value=return_value)
506            ENDIF
507         ENDIF
508         CALL  part_io%sm_node_barrier()
509         CALL cpu_log( log_point_s(89), 'particle NetCDF output', 'stop' )
510
511#if defined( __parallel )
512         CALL MPI_BARRIER(comm2d, ierr)               !kk This Barrier is necessary, not sure why
513#endif
514      END DO
515
516      CALL deallocate_and_free
517
518!     write Time value
519      IF(myid == 0)   THEN
520         ALLOCATE(my_time(1))
521         bounds_start(1)  = icount
522         value_counts(1)  = 1
523         my_time(1)       = simulated_time
524         IF(unlimited_dimension)   THEN               ! For workaround described in dop finalize
525            time_axis_values(icount) = my_time(1)
526         ELSE
527            CALL  netcdf4_write_variable('parallel', file_id, var_id%time, bounds_start(1:1),      &
528                                          value_counts(1:1), bounds_origin(1:1),                   &
529                                         .TRUE., values_realwp_1d=my_time, return_value=return_value)
530         ENDIF
531         DEALLOCATE(my_time)
532      ENDIF
533
534      RETURN
535   END SUBROUTINE dop_output_tseries
536
537   SUBROUTINE dop_finalize
538      IMPLICIT NONE
539#if defined( __netcdf4 )
540      INTEGER(iwp) :: var_len
541#endif
542      INTEGER(iwp) :: return_value                !< Return value data_output_netcdf4 .. routines
543#if defined( __netcdf4 )
544      INTEGER(iwp) :: ierr                        !< MPI error code
545#endif
546
547      IF( win_prt_i /= -1)   THEN
548         CALL part_io%sm_free_shared (win_prt_i)
549      ENDIF
550      IF( win_prt_r /= -1)   THEN
551         CALL part_io%sm_free_shared (win_prt_r)
552      ENDIF
553
554      IF( file_id /= -1 .AND. part_io%iam_io_pe)   THEN
555         CALL netcdf4_finalize( 'parallel', file_id, return_value )
556         file_id = -1
557
558#if defined( __netcdf4 )
559!
560!--      For yet unknown reasons it is not possible to write in parallel mode the time values to NetCDF variable time
561!--      This workaround closes the parallel file and writes the time values in sequential mode on PE0
562!kk      This is a real quick and dirty workaround and the problem should be solved in the final version !!!
563
564         If(myid == 0 .AND. unlimited_dimension)   THEN
565            ierr = nf90_open (TRIM(file_name),NF90_WRITE, file_id)
566            ierr = nf90_inquire_dimension(file_id, did%time, len = var_len)
567
568            ierr = nf90_put_var (file_id, var_id%time, time_axis_values(1:var_len))
569
570            ierr = nf90_close (file_id)
571         ENDIF
572#endif
573      ENDIF
574
575      RETURN
576   END SUBROUTINE dop_finalize
577
578!  Private subroutines
579
580
581!
582!  kk    Not sure if necessary, but set ALL particle numnber to -1
583
584   SUBROUTINE set_indef_particle_nr
585
586      IMPLICIT NONE
587      INTEGER(iwp) :: i                           !<
588      INTEGER(iwp) :: j                           !<
589      INTEGER(iwp) :: k                           !<
590      INTEGER(iwp) :: n                           !<
591
592      DO i=nxl,nxr
593         DO j=nys,nyn
594            DO k=nzb+1,nzt
595               DO n=1,SIZE(grid_particles(k,j,i)%particles)
596                  grid_particles(k,j,i)%particles(n)%particle_nr = -1
597               END DO
598            END DO
599         ENDDO
600      ENDDO
601
602      RETURN
603   END SUBROUTINE set_indef_particle_nr
604!
605!-- Count particles scheduled for output
606!-- here are pts_increment and pts_percentage are used to select output particles
607
608   SUBROUTINE count_output_particles (pcount)
609      IMPLICIT NONE
610
611      INTEGER(iwp), INTENT(OUT)  :: pcount        !<
612
613      INTEGER(iwp) :: i                           !<
614      INTEGER(iwp) :: j                           !<
615      INTEGER(iwp) :: k                           !<
616      INTEGER(iwp) :: n                           !<
617      INTEGER(iwp) :: n_all                       !< count all particles for MOD function
618      REAL(dp)     :: fcount
619      REAL(dp)     :: finc
620
621      pcount = 0
622      IF(pts_increment == 1)   THEN
623         pcount = SUM(prt_count)
624      ELSE
625         n_all = 0
626         DO i=nxl,nxr
627            DO j=nys,nyn
628               DO k=nzb+1,nzt
629                  DO n=1,prt_count(k,j,i)
630                     IF(MOD(n_all,pts_increment) == 0)   THEN
631                        pcount = pcount+1
632                     ENDIF
633                     n_all = n_all+1
634                  END DO
635               END DO
636            ENDDO
637         ENDDO
638      ENDIF
639
640      IF(pts_percentage < 100. )   THEN
641
642         finc   = pts_percentage/100
643
644         pcount = 0
645         fcount = 0.0
646
647         DO i=nxl,nxr
648            DO j=nys,nyn
649               DO k=nzb+1,nzt
650                  DO n=1,prt_count(k,j,i)
651                     fcount = fcount + finc
652                     IF(pcount < int(fcount) )  THEN
653                        pcount = pcount+1
654                     ENDIF
655                  END DO
656               END DO
657            ENDDO
658         ENDDO
659
660      ENDIF
661
662      RETURN
663   END SUBROUTINE count_output_particles
664
665   SUBROUTINE dop_read_output_particle_list (nr_particles_out)
666      IMPLICIT NONE
667      INTEGER(iwp),INTENT(OUT)           :: nr_particles_out
668
669      INTEGER(iwp)             :: i      !<
670      INTEGER(iwp)             :: iu     !<
671      INTEGER(iwp)             :: istat  !<
672      INTEGER(idp)             :: dummy  !<
673
674      iu    = 345
675      istat = 0
676      nr_particles_out = 0
677
678      OPEN(unit=iu, file=TRIM(pts_id_file))      !kk should be changed to check_open
679
680!     First stridem cout output particle
681
682      DO WHILE (istat == 0)
683         READ(iu,*,iostat=istat)  dummy
684         nr_particles_out = nr_particles_out+1
685      END DO
686
687      nr_particles_out = nr_particles_out-1             ! subtract 1 for end of file read
688
689      ALLOCATE(part_id_list_file(nr_particles_out))
690
691      REWIND(iu)
692
693!--   second stride, read particle ids for scheduled output particle
694
695      DO i=1,nr_particles_out
696         READ(iu,*) part_id_list_file(i)
697      END DO
698
699
700      CLOSE (iu)
701
702   END SUBROUTINE dop_read_output_particle_list
703!
704!-- Setb output particle number for selected active particles
705
706   SUBROUTINE set_particle_number
707      IMPLICIT NONE
708
709      INTEGER(iwp) :: i                           !<
710      INTEGER(iwp) :: j                           !<
711      INTEGER(iwp) :: k                           !<
712      INTEGER(iwp) :: n                           !<
713      INTEGER(iwp) :: n_all                       !< count all particles for MOD function
714      INTEGER(iwp) :: particle_nr                 !< output particle number
715      INTEGER(iwp) :: pcount                      !< local particle count in case of pts_percentage
716      REAL(dp)     :: fcount                      !< partical progress in %/100
717      REAL(dp)     :: finc                        !< increment of particle
718
719      pcount = 0
720      fcount = 0.0
721
722      particle_nr = start_local_numbering
723      n_all       = 0
724      DO i=nxl,nxr
725         DO j=nys,nyn
726            DO k=nzb+1,nzt
727               IF(pts_increment > 1)   THEN
728                  DO n=1,prt_count(k,j,i)
729                     IF(MOD(n_all,pts_increment) == 0)   THEN
730                        grid_particles(k,j,i)%particles(n)%particle_nr = particle_nr
731                        particle_nr = particle_nr+1
732                     ELSE
733                        grid_particles(k,j,i)%particles(n)%particle_nr = -2
734                     ENDIF
735                     n_all = n_all+1
736                  END DO
737               ELSE IF(pts_percentage < 100. )   THEN
738                  finc   = pts_percentage/100
739
740                  DO n=1,prt_count(k,j,i)
741!
742!--                  Every particle move fraction on particle axis (i.e part percent == 80; move 0.8)
743!--                  if increases next whole number, the particle is taken for output
744                     fcount = fcount + finc
745                     IF(pcount < int(fcount) )  THEN
746                        pcount = pcount+1
747                        grid_particles(k,j,i)%particles(n)%particle_nr = particle_nr
748                        particle_nr = particle_nr+1
749                     ELSE
750                        grid_particles(k,j,i)%particles(n)%particle_nr = -2
751                     ENDIF
752                  END DO
753               ELSE
754                  DO n=1,prt_count(k,j,i)
755                     grid_particles(k,j,i)%particles(n)%particle_nr = particle_nr
756                     particle_nr = particle_nr+1
757                  END DO
758               ENDIF
759            END DO
760         ENDDO
761      ENDDO
762
763      RETURN
764   END SUBROUTINE set_particle_number
765
766   SUBROUTINE dop_find_particle_in_outlist
767      IMPLICIT NONE
768
769      INTEGER(iwp) :: i                           !<
770#if defined( __parallel )
771      INTEGER(iwp) :: ierr                        !< MPI error code
772#endif
773      INTEGER(iwp) :: j                           !<
774      INTEGER(iwp) :: k                           !<
775      INTEGER(iwp) :: l                           !<
776      INTEGER(iwp) :: n                           !<
777      INTEGER(iwp) :: nr_part                     !<
778!      INTEGER, save :: icount=0
779
780      nr_part = 0
781
782!
783!--   If there is a long particle output list, for performance reason it may become necessary to optimize the
784!--   following loop.
785!
786!--       Decode the particle id in i,j,k,n
787!--       Split serach, i.e search first in i, then in j ...
788
789      DO i=nxl,nxr
790         DO j=nys,nyn
791            DO k=nzb+1,nzt
792               DO n=1,prt_count(k,j,i)
793                  DO l=1,SIZE(part_id_list_file)
794                     IF(grid_particles(k,j,i)%particles(n)%id == part_id_list_file(l))   THEN
795                        grid_particles(k,j,i)%particles(n)%particle_nr =  l
796                        nr_part = nr_part+1
797                     ENDIF
798                  END DO
799               END DO
800            END DO
801         ENDDO
802      ENDDO
803
804#if defined( __parallel )
805      CALL MPI_ALLREDUCE( nr_part, initial_number_of_active_particles, 1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
806#else
807      initial_number_of_active_particles = nr_part
808#endif
809
810   END SUBROUTINE dop_find_particle_in_outlist
811
812!-- Netcdf Setup
813!
814!--   Open NetCDF File DATA_1D_PTS_NETCDF
815!--   Define Dimensions and variables
816!--   Write constant variables
817
818   SUBROUTINE dop_netcdf_setup
819      IMPLICIT NONE
820
821      INTEGER,PARAMETER        :: global_id_in_file = -1
822      INTEGER(iwp)             :: i
823      INTEGER(iwp)             :: fix_ind
824      INTEGER(iwp)             :: var_ind
825
826      INTEGER(iwp)             :: return_value
827      LOGICAL                  :: const_flag
828
829      INTEGER, DIMENSION(2)    :: dimension_ids
830
831      nr_time_values = end_time/dt_dopts+1               !kk has to be adapted to formular of 3d output
832
833      IF (part_io%iam_io_pe)   THEN
834         CALL netcdf4_init_module( "", part_io%comm_io, 0, 9, .TRUE., -1 )
835
836         file_name = 'DATA_1D_PTS_NETCDF'
837#if defined( __parallel )
838         CALL netcdf4_open_file( 'parallel', trim(file_name), file_id, return_value )
839#else
840         CALL netcdf4_open_file( 'serial', trim(file_name), file_id, return_value )
841#endif
842!
843!--      global attributes
844
845         CALL netcdf4_write_attribute( 'parallel', file_id, global_id_in_file, 'comment', &
846                     'Particle ouput created by PALM module data_output_particle', return_value=return_value )
847
848         CALL netcdf4_write_attribute( 'parallel', file_id, global_id_in_file, 'initial_nr_particles', &
849                     value_int32=initial_number_of_active_particles, return_value=return_value )
850!
851!--      define dimensions
852         CALL netcdf4_init_dimension( 'parallel', file_id, did%prt, var_id%prt, &
853               'prt','int32' , nr_particles_out, .TRUE., return_value )
854
855         IF(unlimited_dimension)   THEN
856            CALL netcdf4_init_dimension( 'parallel', file_id, did%time, var_id%time, &
857                                 'time','real32' , -1, .TRUE., return_value )
858            ALLOCATE (time_axis_values(nr_time_values))
859         ELSE
860            CALL netcdf4_init_dimension( 'parallel', file_id, did%time, var_id%time, &
861                                 'time','real32' , nr_time_values, .TRUE., return_value )
862         ENDIF
863      END IF
864!
865!--   Variables without time axis
866!--   These variables will always be written only once at the beginning of the file
867
868      dimension_ids(1) = did%prt
869
870      fix_ind = 1
871      fix_variables(fix_ind)%name  = 'origin_x'
872      fix_variables(fix_ind)%units = 'meter'
873
874      IF (part_io%iam_io_pe)   THEN
875         CALL netcdf4_init_variable( 'parallel', file_id, fix_variables(fix_ind)%var_id, fix_variables(fix_ind)%name, 'real32',  &
876                                   dimension_ids(1:1), .FALSE., return_value )
877
878         CALL netcdf4_write_attribute( 'parallel', file_id, fix_variables(fix_ind)%var_id, 'units', &
879                                   value_char=trim(fix_variables(fix_ind)%units), return_value=return_value )
880      ENDIF
881
882      fix_ind = fix_ind+1
883      fix_variables(fix_ind)%name  = 'origin_y'
884      fix_variables(fix_ind)%units = 'meter'
885
886      IF (part_io%iam_io_pe)   THEN
887         CALL netcdf4_init_variable( 'parallel', file_id, fix_variables(fix_ind)%var_id, fix_variables(fix_ind)%name, 'real32',  &
888                                   dimension_ids(1:1), .FALSE., return_value )
889
890         CALL netcdf4_write_attribute( 'parallel', file_id, fix_variables(fix_ind)%var_id, 'units', &
891                                   value_char=trim(fix_variables(fix_ind)%units), return_value=return_value )
892
893      ENDIF
894
895      fix_ind = fix_ind+1
896      fix_variables(fix_ind)%name  = 'origin_z'
897      fix_variables(fix_ind)%units = 'meter'
898
899      IF (part_io%iam_io_pe)   THEN
900         CALL netcdf4_init_variable( 'parallel', file_id, fix_variables(fix_ind)%var_id, fix_variables(fix_ind)%name, 'real32',  &
901                                   dimension_ids(1:1), .FALSE., return_value )
902
903         CALL netcdf4_write_attribute( 'parallel', file_id, fix_variables(fix_ind)%var_id, 'units', &
904                                   value_char=trim(fix_variables(fix_ind)%units), return_value=return_value )
905      ENDIF
906
907!
908!--   These variables are written if name end with _const'
909      DO i=1,size(data_output_pts)
910         const_flag = (INDEX(TRIM(data_output_pts(i)),'_const') > 0)
911         IF(LEN(TRIM(data_output_pts(i))) > 0 .AND. const_flag)    THEN
912            fix_ind = fix_ind+1
913            fix_variables(fix_ind)%name  = TRIM(data_output_pts(i))
914
915            SELECT CASE (TRIM(fix_variables(fix_ind)%name))
916               CASE('radius_const')
917                  fix_variables(fix_ind)%units = 'meter'
918                  fix_variables(fix_ind)%is_integer = .FALSE.
919               CASE('aux1_const')
920                  fix_variables(fix_ind)%units = 'depend_on_setup'
921                  fix_variables(fix_ind)%is_integer = .FALSE.
922               CASE('aux2_const')
923                  fix_variables(fix_ind)%units = 'depend_on_setup'
924                  fix_variables(fix_ind)%is_integer = .FALSE.
925               CASE('rvar1_const')
926                  fix_variables(fix_ind)%units = 'depend_on_setup'
927                  fix_variables(fix_ind)%is_integer = .FALSE.
928               CASE('rvar2_const')
929                  fix_variables(fix_ind)%units = 'depend_on_setup'
930                  fix_variables(fix_ind)%is_integer = .FALSE.
931               CASE('rvar3_const')
932                  fix_variables(fix_ind)%units = 'depend_on_setup'
933                  fix_variables(fix_ind)%is_integer = .FALSE.
934            END SELECT
935
936            IF (part_io%iam_io_pe)   THEN
937               IF(fix_variables(fix_ind)%is_integer)   THEN
938                  CALL netcdf4_init_variable( 'parallel', file_id, fix_variables(fix_ind)%var_id,  &
939                                              fix_variables(fix_ind)%name, 'int32',                &
940                                              dimension_ids(1:1), .FALSE., return_value )
941               ELSE
942                  CALL netcdf4_init_variable( 'parallel', file_id, fix_variables(fix_ind)%var_id,  &
943                                               fix_variables(fix_ind)%name, 'real32',              &
944                                               dimension_ids(1:1), .FALSE., return_value )
945               ENDIF
946
947               CALL netcdf4_write_attribute( 'parallel', file_id, fix_variables(fix_ind)%var_id, 'units', &
948                  value_char=trim(fix_variables(fix_ind)%units), return_value=return_value )
949            ENDIF
950
951         ENDIF
952      ENDDO
953
954      nr_fix_variables = fix_ind
955!
956!--   Variables time axis
957!--   These variables will always be written in the time loop
958
959      dimension_ids(1) = did%prt
960      dimension_ids(2) = did%time
961
962      var_ind = 0
963
964      DO i=1,size(data_output_pts)
965         const_flag = (INDEX(TRIM(data_output_pts(i)),'_const') > 0)
966         IF(LEN(TRIM(data_output_pts(i))) > 0 .AND. .NOT. const_flag)    THEN
967            var_ind = var_ind+1
968            variables(var_ind)%name  = TRIM(data_output_pts(i))
969
970            SELECT CASE (TRIM(variables(var_ind)%name))
971               CASE('id')
972                  variables(var_ind)%name  = TRIM(data_output_pts(i)) // '_low'
973                  variables(var_ind)%units = 'Number'
974                  variables(var_ind)%is_integer = .TRUE.
975                  IF (part_io%iam_io_pe)   THEN
976                     CALL netcdf4_init_variable( 'parallel', file_id, variables(var_ind)%var_id, variables(var_ind)%name,  &
977                                                'int32', dimension_ids(1:2), .FALSE., return_value )
978                     CALL netcdf4_write_attribute( 'parallel', file_id, variables(var_ind)%var_id, 'units', &
979                                                value_char=trim(variables(var_ind)%units), return_value=return_value )
980                  ENDIF
981
982                  var_ind = var_ind+1
983                  variables(var_ind)%name  = TRIM(data_output_pts(i)) // '_high'
984                  variables(var_ind)%units = 'Number'
985                  variables(var_ind)%is_integer = .TRUE.
986               CASE('particle_nr')
987                  variables(var_ind)%units = 'Number'
988                  variables(var_ind)%is_integer = .TRUE.
989               CASE('class')
990                  variables(var_ind)%units = 'Number'
991                  variables(var_ind)%is_integer = .TRUE.
992               CASE('group')
993                  variables(var_ind)%units = 'Number'
994                  variables(var_ind)%is_integer = .TRUE.
995               CASE('x')
996                  variables(var_ind)%units = 'meter'
997                  variables(var_ind)%is_integer = .FALSE.
998               CASE('y')
999                  variables(var_ind)%units = 'meter'
1000                  variables(var_ind)%is_integer = .FALSE.
1001               CASE('z')
1002                  variables(var_ind)%units = 'meter'
1003                  variables(var_ind)%is_integer = .FALSE.
1004               CASE('speed_x')
1005                  variables(var_ind)%units = 'm/s'
1006                  variables(var_ind)%is_integer = .FALSE.
1007               CASE('speed_y')
1008                  variables(var_ind)%units = 'm/s'
1009                  variables(var_ind)%is_integer = .FALSE.
1010               CASE('speed_z')
1011                  variables(var_ind)%units = 'm/s'
1012                  variables(var_ind)%is_integer = .FALSE.
1013               CASE('radius')
1014                  variables(var_ind)%units = 'meter'
1015                  variables(var_ind)%is_integer = .FALSE.
1016               CASE('age')
1017                  variables(var_ind)%units = 'sec'
1018                  variables(var_ind)%is_integer = .FALSE.
1019               CASE('age_m')
1020                  variables(var_ind)%units = 'sec'
1021                  variables(var_ind)%is_integer = .FALSE.
1022               CASE('dt_sum')
1023                  variables(var_ind)%units = 'sec'
1024                  variables(var_ind)%is_integer = .FALSE.
1025               CASE('e_m')
1026                  variables(var_ind)%units = 'Ws'
1027                  variables(var_ind)%is_integer = .FALSE.
1028               CASE('weight_factor')
1029                  variables(var_ind)%units = 'factor'
1030                  variables(var_ind)%is_integer = .FALSE.
1031               CASE('aux1')
1032                  variables(var_ind)%units = 'depend_on_setup'
1033                  variables(var_ind)%is_integer = .FALSE.
1034               CASE('aux2')
1035                  variables(var_ind)%units = 'depend_on_setup'
1036                  variables(var_ind)%is_integer = .FALSE.
1037               CASE('rvar1')
1038                  variables(var_ind)%units = 'depend_on_setup'
1039                  variables(var_ind)%is_integer = .FALSE.
1040               CASE('rvar2')
1041                  variables(var_ind)%units = 'depend_on_setup'
1042                  variables(var_ind)%is_integer = .FALSE.
1043               CASE('rvar3')
1044                  variables(var_ind)%units = 'depend_on_setup'
1045                  variables(var_ind)%is_integer = .FALSE.
1046            END SELECT
1047
1048            IF (part_io%iam_io_pe)   THEN
1049               IF(variables(var_ind)%is_integer)   THEN
1050                  CALL netcdf4_init_variable( 'parallel', file_id, variables(var_ind)%var_id, variables(var_ind)%name, 'int32',  &
1051                                                  dimension_ids(1:2), .FALSE., return_value )
1052               ELSE
1053                  CALL netcdf4_init_variable( 'parallel', file_id, variables(var_ind)%var_id, variables(var_ind)%name, 'real32',  &
1054                                                  dimension_ids(1:2), .FALSE., return_value )
1055               ENDIF
1056
1057               CALL netcdf4_write_attribute( 'parallel', file_id, variables(var_ind)%var_id, 'units', &
1058                  value_char=trim(variables(var_ind)%units), return_value=return_value )
1059            ENDIF
1060
1061         ENDIF
1062      ENDDO
1063
1064      nr_variables = var_ind
1065
1066      IF (part_io%iam_io_pe)   THEN
1067         CALL netcdf4_stop_file_header_definition( 'parallel', file_id, return_value )
1068      ENDIF
1069
1070      CALL dop_write_axis
1071
1072      RETURN
1073
1074    CONTAINS
1075       SUBROUTINE dop_write_axis
1076          IMPLICIT NONE
1077          INTEGER(iwp)         :: i
1078          INTEGER,DIMENSION(1) :: bounds_start, value_counts, bounds_origin
1079          INTEGER(iwp)         :: return_value                !< Return value data_output_netcdf4 .. routines
1080
1081          INTEGER, POINTER, CONTIGUOUS, DIMENSION(:)   :: prt_val
1082
1083          bounds_origin = 1
1084          bounds_start(1) = 1
1085
1086          IF(myid == 0)   THEN
1087
1088             ALLOCATE(prt_val(nr_particles_out))
1089             DO i=1,nr_particles_out
1090                prt_val(i) = i
1091             ENDDO
1092             value_counts(1) = nr_particles_out
1093
1094             CALL  netcdf4_write_variable('parallel', file_id, var_id%prt, bounds_start, value_counts, bounds_origin, &
1095                .TRUE., values_int32_1d=prt_val, return_value=return_value)
1096
1097             DEALLOCATE(prt_val)
1098          END IF
1099
1100          RETURN
1101       END SUBROUTINE dop_write_axis
1102
1103   END SUBROUTINE dop_netcdf_setup
1104!
1105!- write constant variables
1106
1107   SUBROUTINE dop_write_fixed_variables
1108      IMPLICIT NONE
1109
1110      INTEGER(iwp)             :: i                              !
1111      INTEGER(iwp)             :: return_value                   !
1112
1113      INTEGER(iwp),DIMENSION(1)   :: bounds_origin, bounds_start, value_counts
1114
1115      bounds_origin(1) = 1
1116      bounds_start(1)  = io_start_index
1117      value_counts(1)  = io_end_index-io_start_index+1
1118
1119
1120      DO i=1,nr_fix_variables
1121
1122#if defined( __netcdf4 )
1123         IF(fix_variables(i)%is_integer)   THEN
1124            out_buf_i = NF90_FILL_INT
1125         ELSE
1126            out_buf_r = NF90_FILL_REAL
1127         ENDIF
1128#endif
1129
1130         CALL dop_fill_out_buf (fix_variables(i))
1131
1132         CALL  dop_get_remote_particle (fix_variables(i)%is_integer)
1133
1134         CALL  part_io%sm_node_barrier()
1135         IF (part_io%iam_io_pe)   THEN
1136            IF(fix_variables(i)%is_integer)   THEN
1137               CALL  netcdf4_write_variable('parallel', file_id, fix_variables(i)%var_id, bounds_start, value_counts, &
1138                                            bounds_origin, .FALSE., values_int32_1d=out_buf_i, return_value=return_value)
1139            ELSE
1140               CALL  netcdf4_write_variable('parallel', file_id, fix_variables(i)%var_id, bounds_start, value_counts, &
1141                                            bounds_origin, .FALSE., values_real32_1d=out_buf_r, return_value=return_value)
1142            ENDIF
1143         ENDIF
1144         CALL  part_io%sm_node_barrier()
1145      ENDDO
1146
1147      RETURN
1148   END SUBROUTINE dop_write_fixed_variables
1149
1150   SUBROUTINE dop_delete_particle_number
1151      IMPLICIT NONE
1152
1153      INTEGER(iwp) :: i                           !<
1154      INTEGER(iwp) :: j                           !<
1155      INTEGER(iwp) :: k                           !<
1156      INTEGER(iwp) :: n                           !<
1157
1158!
1159!--   delete inactive particles, i.e. all particles with particle_mask == .FALSE.
1160!--   get an output particle nr of -1
1161!kk   Not sure if it is required here or already done in lagrangian_particle_model_mod before this call
1162
1163      remote_nr_particles = 0
1164      DO i=nxl,nxr
1165         DO j=nys,nyn
1166            DO k=nzb+1,nzt
1167               DO n=1,prt_count(k,j,i)
1168                  IF( .NOT. grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1169                     grid_particles(k,j,i)%particles(n)%particle_nr = -1
1170                  END IF
1171               END DO
1172            END DO
1173         ENDDO
1174      ENDDO
1175
1176      RETURN
1177   END SUBROUTINE dop_delete_particle_number
1178
1179   SUBROUTINE dop_newly_generated_particles
1180      IMPLICIT NONE
1181
1182      INTEGER(iwp) :: i                           !<
1183      INTEGER(iwp) :: j                           !<
1184      INTEGER(iwp) :: k                           !<
1185      INTEGER(iwp) :: n                           !<
1186      INTEGER(iwp) :: n_all                       !< count all particles for MOD function
1187      INTEGER(iwp) :: nr_new_particle             !<
1188      INTEGER(iwp) :: particle_nr                 !<
1189#if defined( __parallel )
1190      INTEGER(iwp) :: ierr                        !< MPI error code
1191#endif
1192      REAL(dp)     :: fcount
1193      REAL(dp)     :: finc
1194      INTEGER(iwp),DIMENSION(0:numprocs-1)     :: nr_particles_new_s
1195      INTEGER(iwp),DIMENSION(0:numprocs-1)     :: nr_particles_new_r
1196      INTEGER(iwp)                             :: start_new_numbering
1197
1198
1199!
1200!--   Count Number of Newly Generated particles
1201!--   Condition for a newparticle: particle_mask = .TRUE. and particle nr of -1
1202!
1203!--   For performance reasons, this subroutine may be combined later with dop_delete_particle_number
1204
1205      nr_new_particle     = 0
1206
1207      IF(pts_increment > 1)   THEN
1208         n_all = 0
1209         DO i=nxl,nxr
1210            DO j=nys,nyn
1211               DO k=nzb+1,nzt
1212                  DO n=1,prt_count(k,j,i)
1213                     IF( grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1214                        IF(grid_particles(k,j,i)%particles(n)%particle_nr == -1)   THEN
1215                           IF(MOD(n_all,pts_increment) == 0)   THEN
1216                              nr_new_particle = nr_new_particle+1
1217                           ENDIF
1218                           n_all = n_all+1
1219                        ENDIF
1220                     END IF
1221                  END DO
1222               END DO
1223            ENDDO
1224         ENDDO
1225      ELSEIF(pts_percentage < 100. )   THEN
1226         finc   = pts_percentage/100
1227         fcount = 0.0
1228
1229         DO i=nxl,nxr
1230            DO j=nys,nyn
1231               DO k=nzb+1,nzt
1232                  DO n=1,prt_count(k,j,i)
1233                     IF( grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1234                        IF(grid_particles(k,j,i)%particles(n)%particle_nr == -1)   THEN
1235                           fcount = fcount + finc
1236                           IF(nr_new_particle < int(fcount) )  THEN
1237                              nr_new_particle = nr_new_particle+1
1238                           ENDIF
1239                        ENDIF
1240                     END IF
1241                  END DO
1242               END DO
1243            ENDDO
1244         ENDDO
1245
1246      ELSE
1247         DO i=nxl,nxr
1248            DO j=nys,nyn
1249               DO k=nzb+1,nzt
1250                  DO n=1,prt_count(k,j,i)
1251                     IF( grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1252                        IF(grid_particles(k,j,i)%particles(n)%particle_nr == -1)   THEN
1253                           nr_new_particle = nr_new_particle+1
1254                        ENDIF
1255                     END IF
1256                  END DO
1257               END DO
1258            ENDDO
1259         ENDDO
1260      ENDIF
1261!
1262!--   Determine start number of new particles on every thread
1263
1264      nr_particles_new_s = 0
1265      nr_particles_new_s(myid) = nr_new_particle
1266
1267#if defined( __parallel )
1268      CALL MPI_ALLREDUCE( nr_particles_new_s, nr_particles_new_r, SIZE( nr_particles_new_s ), MPI_INTEGER,    &
1269                                              MPI_SUM, comm2d, ierr )
1270#else
1271      nr_particles_new_r = nr_particles_new_s
1272#endif
1273!
1274!--   Abortion if selected particles from new particle set would exceed particle axis of output
1275!--   changed by JS
1276      IF ( ( SUM(nr_particles_new_r) + initial_number_of_active_particles)  >  nr_particles_out ) THEN
1277         RETURN
1278      ENDIF
1279
1280      start_new_numbering = initial_number_of_active_particles+1
1281
1282      IF ( myid > 0)   THEN
1283         DO i=1,numprocs-1
1284            start_new_numbering = start_new_numbering+nr_particles_new_r(i-1)
1285            IF(myid == i)  EXIT
1286         ENDDO
1287      END IF
1288
1289      initial_number_of_active_particles = initial_number_of_active_particles+SUM(nr_particles_new_r)
1290
1291      dop_last_active_particle = initial_number_of_active_particles
1292!
1293!--   Set number of new particles
1294
1295      particle_nr         = start_new_numbering
1296      nr_new_particle     = 0
1297
1298      IF(pts_increment > 1)   THEN
1299         n_all = 0
1300         DO i=nxl,nxr
1301            DO j=nys,nyn
1302               DO k=nzb+1,nzt
1303                  DO n=1,prt_count(k,j,i)
1304                     IF( grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1305                        IF(grid_particles(k,j,i)%particles(n)%particle_nr == -1)   THEN
1306                           IF(MOD(n_all,pts_increment) == 0)   THEN
1307                              grid_particles(k,j,i)%particles(n)%particle_nr = particle_nr
1308                              particle_nr = particle_nr+1
1309                           ELSE
1310                              grid_particles(k,j,i)%particles(n)%particle_nr = -2
1311                           ENDIF
1312                           n_all = n_all+1
1313                        ENDIF
1314                     END IF
1315                  END DO
1316               END DO
1317            ENDDO
1318         ENDDO
1319      ELSEIF(pts_percentage < 100. )   THEN
1320         finc   = pts_percentage/100
1321         fcount = 0.0
1322
1323         DO i=nxl,nxr
1324            DO j=nys,nyn
1325               DO k=nzb+1,nzt
1326                  DO n=1,prt_count(k,j,i)
1327                     IF( grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1328                        IF(grid_particles(k,j,i)%particles(n)%particle_nr == -1)   THEN
1329                           fcount = fcount + finc
1330                           IF(nr_new_particle < int(fcount) )  THEN
1331                              grid_particles(k,j,i)%particles(n)%particle_nr = particle_nr
1332                              particle_nr = particle_nr+1
1333                              nr_new_particle = nr_new_particle+1
1334                           ELSE
1335                              grid_particles(k,j,i)%particles(n)%particle_nr = -2
1336                           ENDIF
1337                        ENDIF
1338                     END IF
1339                  END DO
1340               END DO
1341            ENDDO
1342         ENDDO
1343      ELSE
1344         DO i=nxl,nxr
1345            DO j=nys,nyn
1346               DO k=nzb+1,nzt
1347                  DO n=1,prt_count(k,j,i)
1348                     IF( grid_particles(k,j,i)%particles(n)%particle_mask)   THEN
1349                        IF(grid_particles(k,j,i)%particles(n)%particle_nr == -1)   THEN
1350                            grid_particles(k,j,i)%particles(n)%particle_nr = particle_nr
1351                            particle_nr = particle_nr+1
1352                        ENDIF
1353                     END IF
1354                  END DO
1355               END DO
1356            ENDDO
1357         ENDDO
1358
1359
1360      ENDIF
1361
1362
1363
1364      RETURN
1365   END SUBROUTINE dop_newly_generated_particles
1366
1367   SUBROUTINE dop_count_remote_particles
1368      IMPLICIT NONE
1369
1370#if defined( __parallel )
1371      INTEGER(iwp) :: i                           !<
1372      INTEGER(iwp) :: j                           !<
1373      INTEGER(iwp) :: k                           !<
1374      INTEGER(iwp) :: n                           !<
1375      INTEGER(iwp) :: iop                         !<
1376      INTEGER(iwp) :: particle_nr
1377      INTEGER(iwp) :: pe_nr
1378      INTEGER(iwp) :: win_size
1379      INTEGER(iwp) :: ierr                        !< MPI error code
1380      INTEGER(iwp), DIMENSION(0:numprocs-1) :: part_ind
1381
1382!
1383!--   Count remote particles
1384
1385      remote_nr_particles = 0
1386      DO i=nxl,nxr
1387         DO j=nys,nyn
1388            DO k=nzb+1,nzt
1389               DO n=1,prt_count(k,j,i)
1390                  particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1391                  IF ( particle_nr > 0)   THEN
1392                     DO iop=0,numprocs-1                  !kk this loop has to be optimized
1393!
1394!--                     Although the counting is local PE based, the following if is io processor based
1395!--                     because particles in MPI shared memory do not have to be transfered
1396                        IF(particle_nr < io_start_index .OR. particle_nr > io_end_index)   THEN
1397                           IF(particle_nr >= mo_indices(1,iop) .AND. particle_nr <= mo_indices(2,iop))   THEN
1398                              remote_nr_particles(2,iop) = remote_nr_particles(2,iop)+1
1399                           ENDIF
1400                        ENDIF
1401                     ENDDO
1402                  END IF
1403               END DO
1404            END DO
1405         ENDDO
1406      ENDDO
1407
1408      remote_nr_particles(1,0) = 0
1409      DO i=1,numprocs-1
1410         remote_nr_particles(1,i) = remote_nr_particles(1,i-1) + remote_nr_particles(2,i-1)
1411      END DO
1412
1413      win_size = sum(remote_nr_particles(2,:))
1414      CALL dop_alloc_rma_mem (transfer_buffer_i, win_size, win_rma_buf_i)
1415      CALL dop_alloc_rma_mem (transfer_buffer_r, win_size, win_rma_buf_r)
1416
1417      CALL MPI_ALLTOALL (remote_nr_particles, 2, MPI_INTEGER, rma_particles, 2, MPI_INTEGER, comm2d, ierr)
1418
1419!
1420!--   The particles indices are the same for all output variables during one time step
1421!--   therefore, the indices are transfered only once here
1422
1423      part_ind = remote_nr_particles(1,:)
1424      transfer_buffer_i = -9999
1425
1426      DO i=nxl,nxr
1427         DO j=nys,nyn
1428            DO k=nzb+1,nzt
1429
1430               DO n=1,prt_count(k,j,i)
1431                  particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1432                  IF(particle_nr < io_start_index .OR. particle_nr > io_end_index)   THEN
1433                     IF (particle_nr > 0)   THEN
1434                        pe_nr = find_pe_from_particle_nr (particle_nr)
1435                        transfer_buffer_i (part_ind(pe_nr)) = particle_nr
1436                        part_ind(pe_nr) = part_ind(pe_nr)+1
1437                     ENDIF
1438                  ENDIF
1439               END DO
1440            END DO
1441         ENDDO
1442      ENDDO
1443
1444      CALL MPI_Barrier (MPI_COMM_WORLD, ierr)
1445
1446      CALL dop_get_remote_indices
1447
1448#endif
1449      RETURN
1450   END SUBROUTINE dop_count_remote_particles
1451
1452!
1453!- Fill ouput buffer
1454!
1455!    local variables values are copied into output buffer
1456!               local here means local to shared memory group
1457!    remote variable values are copied into transfer buffer
1458!               this is done by all threads
1459
1460   SUBROUTINE dop_fill_out_buf (var)
1461      IMPLICIT NONE
1462
1463      TYPE(var_def),INTENT(IN)    :: var
1464
1465      INTEGER(iwp) :: i                           !<
1466      INTEGER(iwp) :: j                           !<
1467      INTEGER(iwp) :: k                           !<
1468      INTEGER(iwp) :: n                           !<
1469      INTEGER(idp) :: pval                        !<
1470      INTEGER(iwp) :: particle_nr
1471      INTEGER(iwp) :: pe_nr
1472      INTEGER(iwp) :: local_len
1473      CHARACTER(len=32) :: local_name
1474      INTEGER(iwp), DIMENSION(0:numprocs-1) :: part_ind
1475
1476      part_ind = remote_nr_particles(1,:)
1477      transfer_buffer_i = -9998
1478!
1479!--   filling output buffer is the same for variable name and variable name_const
1480!--   therefore set local_name without
1481
1482      local_len = INDEX(TRIM(var%name),'_const')
1483      IF(local_len == 0)   THEN
1484         local_name = var%name
1485      ELSE
1486         local_name = var%name(1:local_len-1)
1487      END IF
1488!
1489!--   In this subroutine the particles are seperated:
1490!
1491!--   All particles which are located in the share memory area of the respective IO thread are copied into
1492!--   the output buffer. The other output particle are copied into the transfer buffer.
1493
1494      SELECT CASE (TRIM(local_name))
1495         CASE('origin_x')
1496            DO i=nxl,nxr
1497               DO j=nys,nyn
1498                  DO k=nzb+1,nzt
1499                     DO n=1,prt_count(k,j,i)
1500                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1501                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1502                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%origin_x
1503                        ELSE IF (particle_nr > 0)   THEN
1504                           pe_nr = find_pe_from_particle_nr (particle_nr)
1505                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%origin_x
1506                           part_ind(pe_nr) = part_ind(pe_nr)+1
1507                        ENDIF
1508                     END DO
1509                  END DO
1510               ENDDO
1511            ENDDO
1512         CASE('origin_y')
1513            DO i=nxl,nxr
1514               DO j=nys,nyn
1515                  DO k=nzb+1,nzt
1516                     DO n=1,prt_count(k,j,i)
1517                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1518                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1519                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%origin_y
1520                        ELSE IF (particle_nr > 0)   THEN
1521                           pe_nr = find_pe_from_particle_nr (particle_nr)
1522                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%origin_y
1523                           part_ind(pe_nr) = part_ind(pe_nr)+1
1524                        ENDIF
1525                     END DO
1526                  END DO
1527               ENDDO
1528            ENDDO
1529         CASE('origin_z')
1530            DO i=nxl,nxr
1531               DO j=nys,nyn
1532                  DO k=nzb+1,nzt
1533                     DO n=1,prt_count(k,j,i)
1534                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1535                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1536                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%origin_z
1537                        ELSE IF (particle_nr > 0)   THEN
1538                           pe_nr = find_pe_from_particle_nr (particle_nr)
1539                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%origin_z
1540                           part_ind(pe_nr) = part_ind(pe_nr)+1
1541                        ENDIF
1542                     END DO
1543                  END DO
1544               ENDDO
1545            ENDDO
1546         CASE('id_low')
1547            DO i=nxl,nxr
1548               DO j=nys,nyn
1549                  DO k=nzb+1,nzt
1550                     DO n=1,prt_count(k,j,i)
1551                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1552                        pval = IBITS(grid_particles(k,j,i)%particles(n)%id,0,32)
1553                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1554                           out_buf_i(particle_nr) = INT(pval,4)
1555                        ELSE IF (particle_nr > 0)   THEN
1556                           pe_nr = find_pe_from_particle_nr (particle_nr)
1557                           transfer_buffer_i (part_ind(pe_nr)) = INT(pval,4)
1558                           part_ind(pe_nr) = part_ind(pe_nr)+1
1559                        ENDIF
1560                     END DO
1561                  END DO
1562               ENDDO
1563            ENDDO
1564         CASE('id_high')
1565            DO i=nxl,nxr
1566               DO j=nys,nyn
1567                  DO k=nzb+1,nzt
1568                     DO n=1,prt_count(k,j,i)
1569                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1570                        pval = IBITS(grid_particles(k,j,i)%particles(n)%id,32,32)
1571                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1572                           out_buf_i(particle_nr) = INT(pval,4)
1573                        ELSE IF (particle_nr > 0)   THEN
1574                           pe_nr = find_pe_from_particle_nr (particle_nr)
1575                           transfer_buffer_i (part_ind(pe_nr)) = INT(pval,4)
1576                           part_ind(pe_nr) = part_ind(pe_nr)+1
1577                        ENDIF
1578                     END DO
1579                  END DO
1580               ENDDO
1581            ENDDO
1582         CASE('particle_nr')
1583            DO i=nxl,nxr
1584               DO j=nys,nyn
1585                  DO k=nzb+1,nzt
1586                     DO n=1,prt_count(k,j,i)
1587                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1588                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1589                           out_buf_i(particle_nr) = grid_particles(k,j,i)%particles(n)%particle_nr
1590                        ELSE IF (particle_nr > 0)   THEN
1591                           pe_nr = find_pe_from_particle_nr (particle_nr)
1592                           transfer_buffer_i (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%particle_nr
1593                           part_ind(pe_nr) = part_ind(pe_nr)+1
1594                        ENDIF
1595                     END DO
1596                  END DO
1597               ENDDO
1598            ENDDO
1599         CASE('class')
1600            DO i=nxl,nxr
1601               DO j=nys,nyn
1602                  DO k=nzb+1,nzt
1603                     DO n=1,prt_count(k,j,i)
1604                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1605                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1606                           out_buf_i(particle_nr) = grid_particles(k,j,i)%particles(n)%class
1607                        ELSE IF (particle_nr > 0)   THEN
1608                           pe_nr = find_pe_from_particle_nr (particle_nr)
1609                           transfer_buffer_i (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%class
1610                           part_ind(pe_nr) = part_ind(pe_nr)+1
1611                        ENDIF
1612                     END DO
1613                  END DO
1614               ENDDO
1615            ENDDO
1616         CASE('group')
1617            DO i=nxl,nxr
1618               DO j=nys,nyn
1619                  DO k=nzb+1,nzt
1620                     DO n=1,prt_count(k,j,i)
1621                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1622                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1623                           out_buf_i(particle_nr) = grid_particles(k,j,i)%particles(n)%group
1624                        ELSE IF (particle_nr > 0)   THEN
1625                           pe_nr = find_pe_from_particle_nr (particle_nr)
1626                           transfer_buffer_i (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%group
1627                           part_ind(pe_nr) = part_ind(pe_nr)+1
1628                        ENDIF
1629                     END DO
1630                  END DO
1631               ENDDO
1632            ENDDO
1633         CASE('x')
1634            DO i=nxl,nxr
1635               DO j=nys,nyn
1636                  DO k=nzb+1,nzt
1637                     DO n=1,prt_count(k,j,i)
1638                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1639                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1640                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%x
1641                        ELSE IF (particle_nr > 0)   THEN
1642                           pe_nr = find_pe_from_particle_nr (particle_nr)
1643                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%x
1644                           part_ind(pe_nr) = part_ind(pe_nr)+1
1645                        ENDIF
1646                     END DO
1647                  END DO
1648               ENDDO
1649            ENDDO
1650         CASE('y')
1651            DO i=nxl,nxr
1652               DO j=nys,nyn
1653                  DO k=nzb+1,nzt
1654                     DO n=1,prt_count(k,j,i)
1655                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1656                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1657                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%y
1658                        ELSE IF (particle_nr > 0)   THEN
1659                           pe_nr = find_pe_from_particle_nr (particle_nr)
1660                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%y
1661                           part_ind(pe_nr) = part_ind(pe_nr)+1
1662                        ENDIF
1663                     END DO
1664                  END DO
1665               ENDDO
1666            ENDDO
1667         CASE('z')
1668            DO i=nxl,nxr
1669               DO j=nys,nyn
1670                  DO k=nzb+1,nzt
1671                     DO n=1,prt_count(k,j,i)
1672                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1673                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1674                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%z
1675                        ELSE IF (particle_nr > 0)   THEN
1676                           pe_nr = find_pe_from_particle_nr (particle_nr)
1677                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%z
1678                           part_ind(pe_nr) = part_ind(pe_nr)+1
1679                        ENDIF
1680                     END DO
1681                  END DO
1682               ENDDO
1683            ENDDO
1684         CASE('speed_x')
1685            DO i=nxl,nxr
1686               DO j=nys,nyn
1687                  DO k=nzb+1,nzt
1688                     DO n=1,prt_count(k,j,i)
1689                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1690                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1691                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%speed_x
1692                        ELSE IF (particle_nr > 0)   THEN
1693                           pe_nr = find_pe_from_particle_nr (particle_nr)
1694                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%speed_x
1695                           part_ind(pe_nr) = part_ind(pe_nr)+1
1696                        ENDIF
1697                     END DO
1698                  END DO
1699               ENDDO
1700            ENDDO
1701         CASE('speed_y')
1702            DO i=nxl,nxr
1703               DO j=nys,nyn
1704                  DO k=nzb+1,nzt
1705                     DO n=1,prt_count(k,j,i)
1706                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1707                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1708                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%speed_y
1709                        ELSE IF (particle_nr > 0)   THEN
1710                           pe_nr = find_pe_from_particle_nr (particle_nr)
1711                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%speed_y
1712                           part_ind(pe_nr) = part_ind(pe_nr)+1
1713                        ENDIF
1714                     END DO
1715                  END DO
1716               ENDDO
1717            ENDDO
1718         CASE('speed_z')
1719            DO i=nxl,nxr
1720               DO j=nys,nyn
1721                  DO k=nzb+1,nzt
1722                     DO n=1,prt_count(k,j,i)
1723                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1724                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1725                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%speed_z
1726                        ELSE IF (particle_nr > 0)   THEN
1727                           pe_nr = find_pe_from_particle_nr (particle_nr)
1728                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%speed_z
1729                           part_ind(pe_nr) = part_ind(pe_nr)+1
1730                        ENDIF
1731                     END DO
1732                  END DO
1733               ENDDO
1734            ENDDO
1735         CASE('radius')
1736            DO i=nxl,nxr
1737               DO j=nys,nyn
1738                  DO k=nzb+1,nzt
1739                     DO n=1,prt_count(k,j,i)
1740                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1741                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1742                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%radius
1743                        ELSE IF (particle_nr > 0)   THEN
1744                           pe_nr = find_pe_from_particle_nr (particle_nr)
1745                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%radius
1746                           part_ind(pe_nr) = part_ind(pe_nr)+1
1747                        ENDIF
1748                     END DO
1749                  END DO
1750               ENDDO
1751            ENDDO
1752         CASE('age')
1753            DO i=nxl,nxr
1754               DO j=nys,nyn
1755                  DO k=nzb+1,nzt
1756                     DO n=1,prt_count(k,j,i)
1757                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1758                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1759                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%age
1760                        ELSE IF (particle_nr > 0)   THEN
1761                           pe_nr = find_pe_from_particle_nr (particle_nr)
1762                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%age
1763                           part_ind(pe_nr) = part_ind(pe_nr)+1
1764                        ENDIF
1765                     END DO
1766                  END DO
1767               ENDDO
1768            ENDDO
1769         CASE('age_m')
1770            DO i=nxl,nxr
1771               DO j=nys,nyn
1772                  DO k=nzb+1,nzt
1773                     DO n=1,prt_count(k,j,i)
1774                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1775                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1776                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%age_m
1777                        ELSE IF (particle_nr > 0)   THEN
1778                           pe_nr = find_pe_from_particle_nr (particle_nr)
1779                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%age_m
1780                           part_ind(pe_nr) = part_ind(pe_nr)+1
1781                        ENDIF
1782                     END DO
1783                  END DO
1784               ENDDO
1785            ENDDO
1786         CASE('dt_sum')
1787            DO i=nxl,nxr
1788               DO j=nys,nyn
1789                  DO k=nzb+1,nzt
1790                     DO n=1,prt_count(k,j,i)
1791                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1792                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1793                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%dt_sum
1794                        ELSE IF (particle_nr > 0)   THEN
1795                           pe_nr = find_pe_from_particle_nr (particle_nr)
1796                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%dt_sum
1797                           part_ind(pe_nr) = part_ind(pe_nr)+1
1798                        ENDIF
1799                     END DO
1800                  END DO
1801               ENDDO
1802            ENDDO
1803         CASE('e_m')
1804            DO i=nxl,nxr
1805               DO j=nys,nyn
1806                  DO k=nzb+1,nzt
1807                     DO n=1,prt_count(k,j,i)
1808                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1809                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1810                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%e_m
1811                        ELSE IF (particle_nr > 0)   THEN
1812                           pe_nr = find_pe_from_particle_nr (particle_nr)
1813                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%e_m
1814                           part_ind(pe_nr) = part_ind(pe_nr)+1
1815                        ENDIF
1816                     END DO
1817                  END DO
1818               ENDDO
1819            ENDDO
1820         CASE('weight_factor')
1821            DO i=nxl,nxr
1822               DO j=nys,nyn
1823                  DO k=nzb+1,nzt
1824                     DO n=1,prt_count(k,j,i)
1825                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1826                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1827                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%weight_factor
1828                        ELSE IF (particle_nr > 0)   THEN
1829                           pe_nr = find_pe_from_particle_nr (particle_nr)
1830                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%weight_factor
1831                           part_ind(pe_nr) = part_ind(pe_nr)+1
1832                        ENDIF
1833                     END DO
1834                  END DO
1835               ENDDO
1836            ENDDO
1837         CASE('aux1')
1838            DO i=nxl,nxr
1839               DO j=nys,nyn
1840                  DO k=nzb+1,nzt
1841                     DO n=1,prt_count(k,j,i)
1842                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1843                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1844                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%aux1
1845                        ELSE IF (particle_nr > 0)   THEN
1846                           pe_nr = find_pe_from_particle_nr (particle_nr)
1847                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%aux1
1848                           part_ind(pe_nr) = part_ind(pe_nr)+1
1849                        ENDIF
1850                     END DO
1851                  END DO
1852               ENDDO
1853            ENDDO
1854         CASE('aux2')
1855            DO i=nxl,nxr
1856               DO j=nys,nyn
1857                  DO k=nzb+1,nzt
1858                     DO n=1,prt_count(k,j,i)
1859                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1860                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1861                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%aux2
1862                        ELSE IF (particle_nr > 0)   THEN
1863                           pe_nr = find_pe_from_particle_nr (particle_nr)
1864                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%aux2
1865                           part_ind(pe_nr) = part_ind(pe_nr)+1
1866                        ENDIF
1867                     END DO
1868                  END DO
1869               ENDDO
1870            ENDDO
1871         CASE('rvar1')
1872            DO i=nxl,nxr
1873               DO j=nys,nyn
1874                  DO k=nzb+1,nzt
1875                     DO n=1,prt_count(k,j,i)
1876                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1877                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1878                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%rvar1
1879                        ELSE IF (particle_nr > 0)   THEN
1880                           pe_nr = find_pe_from_particle_nr (particle_nr)
1881                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%rvar1
1882                           part_ind(pe_nr) = part_ind(pe_nr)+1
1883                        ENDIF
1884                     END DO
1885                  END DO
1886               ENDDO
1887            ENDDO
1888         CASE('rvar2')
1889            DO i=nxl,nxr
1890               DO j=nys,nyn
1891                  DO k=nzb+1,nzt
1892                     DO n=1,prt_count(k,j,i)
1893                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1894                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1895                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%rvar2
1896                        ELSE IF (particle_nr > 0)   THEN
1897                           pe_nr = find_pe_from_particle_nr (particle_nr)
1898                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%rvar2
1899                           part_ind(pe_nr) = part_ind(pe_nr)+1
1900                        ENDIF
1901                     END DO
1902                  END DO
1903               ENDDO
1904            ENDDO
1905         CASE('rvar3')
1906            DO i=nxl,nxr
1907               DO j=nys,nyn
1908                  DO k=nzb+1,nzt
1909                     DO n=1,prt_count(k,j,i)
1910                        particle_nr = grid_particles(k,j,i)%particles(n)%particle_nr
1911                        IF(particle_nr >= io_start_index .AND. particle_nr <= io_end_index)   THEN
1912                           out_buf_r(particle_nr) = grid_particles(k,j,i)%particles(n)%rvar3
1913                        ELSE IF (particle_nr > 0)   THEN
1914                           pe_nr = find_pe_from_particle_nr (particle_nr)
1915                           transfer_buffer_r (part_ind(pe_nr)) = grid_particles(k,j,i)%particles(n)%rvar3
1916                           part_ind(pe_nr) = part_ind(pe_nr)+1
1917                        ENDIF
1918                     END DO
1919                  END DO
1920               ENDDO
1921            ENDDO
1922      END SELECT
1923
1924      RETURN
1925   END SUBROUTINE dop_fill_out_buf
1926
1927#if defined( __parallel )
1928!
1929!- Get indices (displacement) of remot particles
1930   SUBROUTINE dop_get_remote_indices
1931
1932      IMPLICIT NONE
1933
1934      INTEGER(iwp) :: i                           !<
1935
1936      INTEGER(iwp) :: bufsize                     !< size of remote indices array
1937      INTEGER(iwp) :: ind_local                   !< index in remore indices array
1938      INTEGER(iwp) :: ierr                        !< MPI error code
1939      INTEGER(KIND=MPI_ADDRESS_KIND)   :: disp    !< displacement in RMA window
1940
1941      bufsize = SUM(rma_particles(2,:))
1942      ALLOCATE (remote_indices(0:bufsize))
1943      remote_indices = -1
1944
1945      ind_local = 0
1946      CALL MPI_WIN_FENCE( 0, win_rma_buf_i, ierr )
1947      DO i=0,numprocs-1
1948         IF (rma_particles(2,i) > 0)   THEN
1949            disp = rma_particles(1,i)
1950            IF(rma_particles(2,i) > 0)   THEN
1951               CALL MPI_GET (remote_indices(ind_local), rma_particles(2,i), MPI_INTEGER, i, disp,     &
1952                                           rma_particles(2,i), MPI_INTEGER,  win_rma_buf_i, ierr)
1953               ind_local = ind_local + rma_particles(2,i)
1954            END IF
1955         END IF
1956      ENDDO
1957      CALL MPI_WIN_FENCE( 0, win_rma_buf_i, ierr )
1958
1959      RETURN
1960   END SUBROUTINE dop_get_remote_indices
1961#endif
1962
1963
1964   SUBROUTINE dop_get_remote_particle (is_integer)
1965
1966      IMPLICIT NONE
1967
1968      LOGICAL,INTENT(IN)   :: is_integer
1969
1970#if defined( __parallel )
1971      INTEGER(iwp) :: i                           !<
1972      INTEGER(iwp) :: j                           !<
1973      INTEGER(iwp) :: bufsize                     !< size of remote data array
1974      INTEGER(iwp) :: ind_local                   !< index in remore indices array
1975      INTEGER(iwp) :: particle_nr                 !< particle number
1976      INTEGER(iwp) :: ierr                        !< MPI error code
1977      INTEGER(KIND=MPI_ADDRESS_KIND)           :: disp        !< displacement in RMA window
1978      REAL(sp),ALLOCATABLE, DIMENSION(:)       :: rma_buf_r   !< buffer to receive remote data (REAL)
1979      INTEGER(iwp),ALLOCATABLE, DIMENSION(:)   :: rma_buf_i   !< buffer to receive remote data (INTEGER)
1980
1981      bufsize   = sum(rma_particles(2,:))
1982      ind_local = 0
1983      ALLOCATE (rma_buf_r(0:bufsize-1))
1984      ALLOCATE (rma_buf_i(0:bufsize-1))
1985
1986      IF(is_integer)   THEN
1987         CALL MPI_WIN_FENCE( 0, win_rma_buf_i, ierr )
1988      ELSE
1989         CALL MPI_WIN_FENCE( 0, win_rma_buf_r, ierr )
1990      ENDIF
1991      DO i=0,numprocs-1
1992         IF (rma_particles(2,i) > 0)   THEN
1993            IF(is_integer)   THEN
1994               disp = rma_particles(1,i)
1995               CALL MPI_GET (rma_buf_i(ind_local), rma_particles(2,i), MPI_INTEGER, i, disp,     &
1996                                        rma_particles(2,i), MPI_INTEGER,  win_rma_buf_i, ierr)
1997               ind_local = ind_local + rma_particles(2,i)
1998            ELSE
1999               disp = rma_particles(1,i)
2000               CALL MPI_GET (rma_buf_r(ind_local), rma_particles(2,i), MPI_real, i, disp,     &
2001                                         rma_particles(2,i), MPI_real,  win_rma_buf_r, ierr)
2002               ind_local = ind_local + rma_particles(2,i)
2003            ENDIF
2004         END IF
2005      ENDDO
2006      IF(is_integer)   THEN
2007         CALL MPI_WIN_FENCE( 0, win_rma_buf_i, ierr )
2008      ELSE
2009         CALL MPI_WIN_FENCE( 0, win_rma_buf_r, ierr )
2010      ENDIF
2011
2012      ind_local = 0
2013
2014      DO i=0,numprocs-1
2015         IF (rma_particles(2,i) > 0)   THEN
2016            IF(is_integer)   THEN
2017!
2018!--            Copy data from remote PEs into output array
2019               DO j=0,rma_particles(2,i)-1
2020                  particle_nr = remote_indices(ind_local)
2021                  out_buf_i(particle_nr) = rma_buf_i(ind_local)
2022                  ind_local = ind_local+1
2023               END DO
2024            ELSE
2025!
2026!--            Copy data from remote PEs into output array
2027
2028               DO j=0,rma_particles(2,i)-1
2029                  particle_nr = remote_indices(ind_local)
2030                  out_buf_r(particle_nr) = rma_buf_r(ind_local)
2031                  ind_local = ind_local+1
2032               END DO
2033            ENDIF
2034         END IF
2035      ENDDO
2036
2037      IF(ALLOCATED(rma_buf_r)) DEALLOCATE(rma_buf_r)
2038      IF(ALLOCATED(rma_buf_i)) DEALLOCATE(rma_buf_i)
2039#else
2040   IF (is_integer)  THEN
2041   ENDIF
2042#endif
2043
2044      RETURN
2045   END SUBROUTINE dop_get_remote_particle
2046
2047#if defined( __parallel )
2048!
2049!- Allocate memory and cread window for one-sided communication (INTEGER 1-D array)
2050   SUBROUTINE dop_alloc_rma_mem_i1( array, idim1, win )
2051      IMPLICIT NONE
2052
2053      INTEGER(isp), DIMENSION(:), POINTER, INTENT(INOUT)   ::  array  !<
2054      INTEGER(iwp), INTENT(IN)                             ::  idim1   !<
2055      INTEGER(iwp), INTENT(OUT)                            ::  win     !<
2056
2057      INTEGER(KIND=MPI_ADDRESS_KIND) ::  winsize  !< size of RMA window
2058      INTEGER                        ::  ierr     !< MPI error code
2059
2060      winsize = max(idim1, 2)
2061
2062      ALLOCATE(array(0:winsize-1))
2063
2064      winsize = winsize * isp
2065
2066      CALL MPI_WIN_CREATE( array, winsize, isp, MPI_INFO_NULL, comm2d, win, ierr )
2067
2068      array = -1
2069
2070      CALL MPI_WIN_FENCE( 0, win, ierr )
2071
2072   END SUBROUTINE dop_alloc_rma_mem_i1
2073#endif
2074
2075#if defined( __parallel )
2076!
2077!- Allocate memory and cread window for one-sided communication (REAL 1-D array)
2078   SUBROUTINE dop_alloc_rma_mem_r1( array, idim1, win )
2079      IMPLICIT NONE
2080
2081      REAL(sp), DIMENSION(:), POINTER, INTENT(INOUT)     ::  array   !<
2082      INTEGER(iwp), INTENT(IN)                           ::  idim1   !<
2083      INTEGER(iwp), INTENT(OUT)                          ::  win     !<
2084
2085      INTEGER(KIND=MPI_ADDRESS_KIND) ::  winsize  !< size of RMA window
2086      INTEGER                        ::  ierr     !< MPI error code
2087
2088
2089      winsize = max(idim1, 2)
2090
2091
2092      ALLOCATE(array(0:winsize-1))
2093
2094      winsize = winsize * sp
2095
2096      CALL MPI_WIN_CREATE( array, winsize, sp, MPI_INFO_NULL, comm2d, win, ierr )
2097
2098      array = -1.0
2099
2100      CALL MPI_WIN_FENCE( 0, win, ierr )
2101
2102   END SUBROUTINE dop_alloc_rma_mem_r1
2103#endif
2104
2105
2106   SUBROUTINE deallocate_and_free
2107      IMPLICIT NONE
2108
2109#if defined( __parallel )
2110      INTEGER                        ::  ierr     !< MPI error code
2111#endif
2112
2113#if defined( __parallel )
2114      CALL MPI_Win_free (win_rma_buf_i, ierr)
2115      CALL MPI_Win_free (win_rma_buf_r, ierr)
2116#endif
2117      IF (ALLOCATED(remote_indices)) DEALLOCATE(remote_indices)
2118
2119      DEALLOCATE(transfer_buffer_i)
2120      DEALLOCATE(transfer_buffer_r)
2121
2122      RETURN
2123
2124   END SUBROUTINE deallocate_and_free
2125
2126   FUNCTION find_pe_from_particle_nr (particle_nr)  RESULT (pe_nr)
2127      IMPLICIT NONE
2128
2129      INTEGER(iwp), INTENT(IN)       :: particle_nr
2130      INTEGER(iwp)                   :: pe_nr             !<
2131      INTEGER(iwp)                   :: base              !<
2132      INTEGER(iwp)                   :: pnr               !<
2133
2134      IF(irregular_distribubtion)   THEN
2135         IF(particle_nr <= nr_particles_rest*nr_particles_PE)   THEN
2136            pe_nr = (particle_nr-1)/nr_particles_PE
2137         ELSE
2138            base  = nr_particles_rest*nr_particles_PE
2139            pnr   = particle_nr - base
2140            pe_nr = (pnr-1)/(nr_particles_PE-1)
2141            pe_nr = pe_nr+nr_particles_rest
2142         ENDIF
2143      ELSE
2144         pe_nr = (particle_nr-1)/nr_particles_PE
2145      ENDIF
2146
2147
2148!kk   This error test is to detect programming errors. For performance reasons it can be removed in
2149!kk   the final, stable version
2150
2151   END FUNCTION find_pe_from_particle_nr
2152
2153END MODULE data_output_particle_mod
Note: See TracBrowser for help on using the repository browser.