Home

Context Navigation

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

Last change on this file since 4785 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

Rev	Line
[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	!
	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 monakurppa $
[4779]	26	! Avoid overlong lines and unused variables
	27	!
	28	! 4778 2020-11-09 13:40:05Z raasch
[4778]	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
[4779]	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)
[4778]	502	ELSE
[4779]	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)
[4778]	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
[4779]	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)
[4778]	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
[4779]	539	#if defined( __netcdf4 )
[4778]	540	INTEGER(iwp) :: var_len
[4779]	541	#endif
[4778]	542	INTEGER(iwp) :: return_value !< Return value data_output_netcdf4 .. routines
[4779]	543	#if defined( __netcdf4 )
[4778]	544	INTEGER(iwp) :: ierr !< MPI error code
[4779]	545	#endif
[4778]	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
[4779]	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 )
[4778]	941	ELSE
[4779]	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 )
[4778]	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
[4779]	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)
[4778]	1139	ELSE
[4779]	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
[4778]	1143	ENDIF
	1144	CALL part_io%sm_node_barrier()
[4779]	1145	ENDDO
[4778]	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 |