Home

Context Navigation

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

Last change on this file since 4783 was 4779, checked in by suehring, 4 years ago
Avoid overlong lines and unused variables
Property svn:keywords set to `Id`
File size: 82.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |