source: palm/trunk/SOURCE/data_output_particle_mod.f90

Last change on this file was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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