Home

Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

restart_data_mpi_io_mod.f90

Last change on this file was 4896, checked in by raasch, 3 years ago
small re-formatting to follow the coding standard, typo in file appendix removed, more meaningful variable names assigned, redundant code removed
Property svn:keywords set to `Id`
File size: 169.1 KB

Line
1	!> @file restart_data_mpi_io_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-2021 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------------------------!
18	!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: restart_data_mpi_io_mod.f90 4896 2021-03-03 16:10:18Z banzhafs $
27	! more meaningful variable names assigned,
28	! redundant code removed
29	!
30	! 4893 2021-03-02 16:39:14Z raasch
31	! revised output of surface data via MPI-IO for better performance
32	!
33	! 4857 2021-01-26 07:24:41Z raasch
34	! bugfix: allocation of 3d-int4 array moved from particle output to standard output
35	!
36	! 4832 2021-01-07 07:15:12Z raasch
37	! some debug output flushed
38	!
39	! 4828 2021-01-05 11:21:41Z Giersch
40	! interface for 3d logical arrays added
41	!
42	! 4694 2020-09-23 15:09:19Z pavelkrc
43	! Fix name of subroutine in message
44	!
45	! 4629 2020-07-29 09:37:56Z raasch
46	! support for MPI Fortran77 interface (mpif.h) removed
47	!
48	! 4628 2020-07-29 07:23:03Z raasch
49	! extensions required for MPI-I/O of particle data to restart files
50	!
51	! 4622 2020-07-23 09:02:23Z raasch
52	! unused variables removed
53	!
54	! 4621 2020-07-23 08:15:59Z raasch
55	! bugfixes for serial (non-parallel) mode
56	!
57	! 4619 2020-07-22 13:21:28Z raasch
58	! unused variable removed
59	!
60	! 4617 2020-07-22 09:48:50Z raasch
61	! Cyclic fill mode implemented
62	!
63	! 4598 2020-07-10 10:13:23Z suehring
64	! Bugfix in treatment of 3D soil arrays
65	!
66	! 4591 2020-07-06 15:56:08Z raasch
67	! File re-formatted to follow the PALM coding standard
68	!
69	! 4539 2020-05-18 14:05:17Z raasch
70	! Checks added, if index limits in header are exceeded
71	! Bugfix in rrd_mpi_io_int_2d
72	!
73	! 4536 2020-05-17 17:24:13Z raasch
74	! Messages and debug output converted to PALM routines
75	!
76	! 4534 2020-05-14 18:35:22Z raasch
77	! I/O on reduced number of cores added (using shared memory MPI)
78	!
79	! 4500 2020-04-17 10:12:45Z suehring
80	! Fix too long lines
81	!
82	! 4498 2020-04-15 14:26:31Z raasch
83	! Bugfix for creation of filetypes, argument removed from rd_mpi_io_open
84	!
85	! 4497 2020-04-15 10:20:51Z raasch
86	! Last bugfix deactivated because of compile problems
87	!
88	! 4496 2020-04-15 08:37:26Z raasch
89	! Problem with posix read arguments for surface data fixed
90	!
91	! 4495 2020-04-13 20:11:20Z raasch
92	! Initial version (K. Ketelsen), adjusted to PALM formatting standards (S. Raasch)
93	!
94	!
95	! Description:
96	! ------------
97	!> Routines for restart data handling using MPI-IO.
98	!--------------------------------------------------------------------------------------------------!
99	MODULE restart_data_mpi_io_mod
100
101	#if defined( __parallel )
102	USE MPI
103	#else
104	USE posix_interface, &
105	ONLY: posix_close, &
106	posix_lseek, &
107	posix_open, &
108	posix_read, &
109	posix_write
110	#endif
111
112	USE, INTRINSIC :: ISO_C_BINDING
113
114	USE control_parameters, &
115	ONLY: debug_output, &
116	debug_string, &
117	include_total_domain_boundaries, &
118	message_string, &
119	restart_data_format_input, &
120	restart_data_format_output, &
121	restart_file_size
122
123	USE exchange_horiz_mod, &
124	ONLY: exchange_horiz, &
125	exchange_horiz_2d
126
127	USE indices, &
128	ONLY: nbgp, &
129	nnx, &
130	nny, &
131	nx, &
132	nxl, &
133	nxlg, &
134	nx_on_file, &
135	nxr, &
136	nxrg, &
137	ny, &
138	nyn, &
139	nyng, &
140	ny_on_file, &
141	nys, &
142	nysg, &
143	nz, &
144	nzb, &
145	nzt
146
147	USE kinds
148
149	USE particle_attributes, &
150	ONLY: grid_particles, &
151	particles, &
152	particle_type, &
153	prt_count, &
154	zero_particle
155
156	USE pegrid, &
157	ONLY: comm1dx, &
158	comm1dy, &
159	comm2d, &
160	communicator_configurations, &
161	myid, &
162	myidx, &
163	myidy, &
164	numprocs, &
165	pdims
166
167	USE shared_memory_io_mod, &
168	ONLY: domain_decomposition_grid_features, &
169	sm_class
170
171
172	IMPLICIT NONE
173
174	CHARACTER(LEN=128) :: io_file_name !> internal variable to communicate filename between different subroutines
175
176	#if defined( __parallel )
177	INTEGER(iwp) :: ierr !< error status of MPI-calls
178	INTEGER(iwp), PARAMETER :: rd_offset_kind = MPI_OFFSET_KIND !< Adress or Offset kind
179	INTEGER(iwp), PARAMETER :: rd_status_size = MPI_STATUS_SIZE !<
180	#else
181	INTEGER(iwp), PARAMETER :: rd_offset_kind = C_SIZE_T !<
182	INTEGER(iwp), PARAMETER :: rd_status_size = 1 !< Not required in sequential mode
183	#endif
184
185	INTEGER(iwp) :: debug_level = 1 !< TODO: replace with standard debug output steering
186
187	INTEGER(iwp) :: comm_io !< Communicator for MPI-IO
188	INTEGER(iwp) :: fh = -1 !< MPI-IO file handle
189	#if defined( __parallel )
190	INTEGER(iwp) :: ft_surf = -1 !< MPI filetype surface data
191	INTEGER(iwp) :: ft_2di_nb !< MPI filetype 2D array INTEGER no outer boundary
192	INTEGER(iwp) :: ft_2d !< MPI filetype 2D array REAL with outer boundaries
193	INTEGER(iwp) :: ft_3d !< MPI filetype 3D array REAL with outer boundaries
194	INTEGER(iwp) :: ft_3di4 = -1 !< MPI filetype 3D array INTEGER*4
195	INTEGER(iwp) :: ft_3di8 = -1 !< MPI filetype 3D array INTEGER*8
196	INTEGER(iwp) :: ft_3dsoil !< MPI filetype for 3d-soil array
197	#endif
198	INTEGER(iwp) :: glo_start !< global start index on this PE
199	#if defined( __parallel )
200	INTEGER(iwp) :: win_2di !<
201	INTEGER(iwp) :: win_2dr !<
202	INTEGER(iwp) :: win_3di4 = -1 !<
203	INTEGER(iwp) :: win_3di8 = -1 !<
204	INTEGER(iwp) :: win_3dr !<
205	INTEGER(iwp) :: win_3ds !<
206	INTEGER(iwp) :: win_end = -1 !<
207	INTEGER(iwp) :: win_glost = -1 !<
208	INTEGER(iwp) :: win_out = -1 !<
209	INTEGER(iwp) :: win_start = -1 !<
210	INTEGER(iwp) :: win_surf = -1 !<
211	#endif
212	INTEGER(iwp) :: total_number_of_surface_elements !< total number of surface elements for one variable
213
214	INTEGER(KIND=rd_offset_kind) :: array_position !<
215	INTEGER(KIND=rd_offset_kind) :: header_position !<
216
217	INTEGER(iwp), DIMENSION(:,:), POINTER, CONTIGUOUS :: array_2di !<
218
219	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: e_end_index !< extended end index, every grid cell has at least one value
220	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: e_start_index !<
221	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: m_end_index !< module copy of end_index
222	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: m_start_index !<
223	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nr_surfaces_in_tb !< number of surface elements in transfer buffer
224	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: s_index_in_tb !< start index of surface elements in transfer buffer
225	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: s_index_in_window !< start index for the pe in the rma window
226	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: transfer_index !<
227	!
228	!-- Indices for cyclic fill
229	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: o_start_index !<
230	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: c_start_index !<
231	!#if defined( __parallel )
232	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: o_end_index !< extended end index, every grid cell has at least one value
233	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: c_end_index !< extended end index, every grid cell has at least one value
234	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: c_global_start !<
235	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: c_global_end !<
236	!#endif
237
238	INTEGER(isp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: array_3di4 !<
239	INTEGER(idp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: array_3di8 !<
240
241	LOGICAL :: filetypes_created !<
242	LOGICAL :: io_on_limited_cores_per_node !< switch to shared memory MPI-IO
243	LOGICAL :: rd_flag !< file is opened for read
244	LOGICAL :: wr_flag !< file is opened for write
245
246	#if defined( __parallel )
247	INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: local_indices
248	#endif
249
250	REAL(wp), DIMENSION(:), POINTER, CONTIGUOUS :: array_out !<
251	#if defined( __parallel )
252	REAL(wp), DIMENSION(:), POINTER, CONTIGUOUS :: array_1d !<
253	#endif
254	REAL(wp), DIMENSION(:,:), POINTER, CONTIGUOUS :: array_2d !<
255	REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: array_3d !<
256	REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS :: array_3d_soil !<
257
258	!
259	!-- Variable to store the grid features (index bounds) of the temporary arrays that are used
260	!-- to read and write the restart data. They differ depending on if the outer boundary of the
261	!-- total domain is contained in the restart data or not. iog stands for IO-grid.
262	TYPE(domain_decomposition_grid_features) :: iog !<
263
264	!
265	!-- General Header (first 32 byte in restart file)
266	TYPE general_header
267	INTEGER(iwp) :: endian !< little endian (1) or big endian (2) internal format
268	INTEGER(iwp) :: i_outer_bound !< if 1, outer boundaries are stored in restart file
269	INTEGER(iwp) :: nr_arrays !< number of arrays in restart files
270	INTEGER(iwp) :: nr_char !< number of Text strings entries in header
271	INTEGER(iwp) :: nr_int !< number of INTEGER entries in header
272	INTEGER(iwp) :: nr_real !< number of REAL entries in header
273	INTEGER(iwp) :: pes_along_x !< number of PEs along x-direction during writing restart file
274	INTEGER(iwp) :: pes_along_y !< number of PEs along y-direction during writing restart file
275	INTEGER(iwp) :: total_nx !< total number of points in x-direction
276	INTEGER(iwp) :: total_ny !< total number of points in y-direction
277	END TYPE general_header
278
279	TYPE(general_header), TARGET, PUBLIC :: tgh !<
280
281	TYPE(sm_class) :: sm_io !<
282
283	!
284	!-- Declaration of varibales for file header section
285	INTEGER(KIND=rd_offset_kind) :: header_int_index
286	INTEGER, PARAMETER :: max_nr_int=256
287	CHARACTER(LEN=32), DIMENSION(max_nr_int) :: int_names
288	INTEGER(KIND=iwp), DIMENSION(max_nr_int) :: int_values
289
290	INTEGER(KIND=rd_offset_kind) :: header_char_index
291	INTEGER, PARAMETER :: max_nr_char=128
292	CHARACTER(LEN=128), DIMENSION(max_nr_char) :: text_lines
293
294	INTEGER(KIND=rd_offset_kind) :: header_real_index
295	INTEGER, PARAMETER :: max_nr_real=256
296	CHARACTER(LEN=32), DIMENSION(max_nr_real) :: real_names
297	REAL(KIND=wp), DIMENSION(max_nr_real) :: real_values
298
299	INTEGER(KIND=rd_offset_kind) :: header_array_index
300	INTEGER, PARAMETER :: max_nr_arrays=600
301	CHARACTER(LEN=32), DIMENSION(max_nr_arrays) :: array_names
302	INTEGER(KIND=rd_offset_kind), DIMENSION(max_nr_arrays) :: array_offset
303
304	!
305	!-- Variables to handle the cyclic fill initialization mode
306	INTEGER :: comm_cyclic_fill !< communicator for cyclic fill PEs
307	#if defined( __parallel )
308	INTEGER :: rmawin_2di !< RMA window 2d INTEGER
309	INTEGER :: rmawin_2d !< RMA window 2d REAL
310	INTEGER :: rmawin_3d !< RMA window 3d
311	#endif
312
313	INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:) :: remote_pe
314	INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:) :: remote_pe_s
315	INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:) :: rma_offset
316	INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:) :: rma_offset_s
317	INTEGER(iwp), ALLOCATABLE, DIMENSION(:,:) :: rmabuf_2di
318
319	LOGICAL :: cyclic_fill_mode !< arrays are filled cyclically with data from prerun
320	LOGICAL :: pe_active_for_read = .TRUE. !< this PE is active for reading data from prerun or
321	!< restart run. For restarts all PEs are active.
322
323	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: rmabuf_2d
324	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: rmabuf_3d
325
326	TYPE(domain_decomposition_grid_features) :: mainrun_grid !< grid variables for the main run
327	TYPE(domain_decomposition_grid_features) :: prerun_grid !< grid variables for the prerun
328
329	!
330	!-- MPI_INTEGER8 is not standard MPI, but is supported on most MPI distibutions
331	!-- If not suppported, a workaround could be enabled with the following preprocessor directive
332	!#if defined( __NO_INTEGER8)
333	! INTEGER :: MPI_INTEGER8 !< MPI data type INTEGER8
334	!#endif
335
336	SAVE
337
338	PRIVATE
339
340	PUBLIC restart_file_size, total_number_of_surface_elements
341
342	!
343	!-- PALM interfaces
344	INTERFACE rd_mpi_io_check_array
345	MODULE PROCEDURE rd_mpi_io_check_array
346	END INTERFACE rd_mpi_io_check_array
347
348	INTERFACE rd_mpi_io_close
349	MODULE PROCEDURE rd_mpi_io_close
350	END INTERFACE rd_mpi_io_close
351
352	INTERFACE rd_mpi_io_check_open
353	MODULE PROCEDURE rd_mpi_io_check_open
354	END INTERFACE rd_mpi_io_check_open
355
356	INTERFACE rd_mpi_io_open
357	MODULE PROCEDURE rd_mpi_io_open
358	END INTERFACE rd_mpi_io_open
359
360	INTERFACE rrd_mpi_io
361	MODULE PROCEDURE rrd_mpi_io_char
362	MODULE PROCEDURE rrd_mpi_io_int
363	MODULE PROCEDURE rrd_mpi_io_int_2d
364	MODULE PROCEDURE rrd_mpi_io_int4_3d
365	MODULE PROCEDURE rrd_mpi_io_int8_3d
366	MODULE PROCEDURE rrd_mpi_io_logical
367	MODULE PROCEDURE rrd_mpi_io_logical_3d
368	MODULE PROCEDURE rrd_mpi_io_real
369	MODULE PROCEDURE rrd_mpi_io_real_2d
370	MODULE PROCEDURE rrd_mpi_io_real_3d
371	MODULE PROCEDURE rrd_mpi_io_real_3d_soil
372	END INTERFACE rrd_mpi_io
373
374	INTERFACE rrd_mpi_io_global_array
375	MODULE PROCEDURE rrd_mpi_io_global_array_int_1d
376	MODULE PROCEDURE rrd_mpi_io_global_array_real_1d
377	MODULE PROCEDURE rrd_mpi_io_global_array_real_2d
378	MODULE PROCEDURE rrd_mpi_io_global_array_real_3d
379	MODULE PROCEDURE rrd_mpi_io_global_array_real_4d
380	END INTERFACE rrd_mpi_io_global_array
381
382	INTERFACE rrd_mpi_io_surface
383	MODULE PROCEDURE rrd_mpi_io_surface
384	MODULE PROCEDURE rrd_mpi_io_surface_2d
385	END INTERFACE rrd_mpi_io_surface
386
387	INTERFACE rrd_mpi_io_particles
388	MODULE PROCEDURE rrd_mpi_io_particles
389	END INTERFACE rrd_mpi_io_particles
390
391	INTERFACE rd_mpi_io_particle_filetypes
392	MODULE PROCEDURE rd_mpi_io_particle_filetypes
393	END INTERFACE rd_mpi_io_particle_filetypes
394
395	INTERFACE rd_mpi_io_surface_filetypes
396	MODULE PROCEDURE rd_mpi_io_surface_filetypes
397	END INTERFACE rd_mpi_io_surface_filetypes
398
399	INTERFACE wrd_mpi_io
400	MODULE PROCEDURE wrd_mpi_io_char
401	MODULE PROCEDURE wrd_mpi_io_int
402	MODULE PROCEDURE wrd_mpi_io_int_2d
403	MODULE PROCEDURE wrd_mpi_io_int4_3d
404	MODULE PROCEDURE wrd_mpi_io_int8_3d
405	MODULE PROCEDURE wrd_mpi_io_logical
406	MODULE PROCEDURE wrd_mpi_io_logical_3d
407	MODULE PROCEDURE wrd_mpi_io_real
408	MODULE PROCEDURE wrd_mpi_io_real_2d
409	MODULE PROCEDURE wrd_mpi_io_real_3d
410	MODULE PROCEDURE wrd_mpi_io_real_3d_soil
411	END INTERFACE wrd_mpi_io
412
413	INTERFACE wrd_mpi_io_global_array
414	MODULE PROCEDURE wrd_mpi_io_global_array_int_1d
415	MODULE PROCEDURE wrd_mpi_io_global_array_real_1d
416	MODULE PROCEDURE wrd_mpi_io_global_array_real_2d
417	MODULE PROCEDURE wrd_mpi_io_global_array_real_3d
418	MODULE PROCEDURE wrd_mpi_io_global_array_real_4d
419	END INTERFACE wrd_mpi_io_global_array
420
421	INTERFACE wrd_mpi_io_particles
422	MODULE PROCEDURE wrd_mpi_io_particles
423	END INTERFACE wrd_mpi_io_particles
424
425	INTERFACE wrd_mpi_io_surface
426	MODULE PROCEDURE wrd_mpi_io_surface
427	MODULE PROCEDURE wrd_mpi_io_surface_2d
428	END INTERFACE wrd_mpi_io_surface
429
430	PUBLIC rd_mpi_io_check_array, &
431	rd_mpi_io_check_open, &
432	rd_mpi_io_close, &
433	rd_mpi_io_open, &
434	rd_mpi_io_particle_filetypes, &
435	rd_mpi_io_surface_filetypes, &
436	rrd_mpi_io, &
437	rrd_mpi_io_global_array, &
438	rrd_mpi_io_particles, &
439	rrd_mpi_io_surface, &
440	wrd_mpi_io, &
441	wrd_mpi_io_global_array, &
442	wrd_mpi_io_particles, &
443	wrd_mpi_io_surface
444
445	CONTAINS
446
447
448	!--------------------------------------------------------------------------------------------------!
449	! Description:
450	! ------------
451	!> Open restart file for read or write with MPI-IO
452	!--------------------------------------------------------------------------------------------------!
453	SUBROUTINE rd_mpi_io_open( action, file_name, open_for_global_io_only )
454
455	IMPLICIT NONE
456
457	CHARACTER(LEN=*), INTENT(IN) :: action !<
458	CHARACTER(LEN=*), INTENT(IN) :: file_name !<
459
460	INTEGER(iwp) :: i !<
461	INTEGER(iwp) :: gh_size !<
462
463	INTEGER(KIND=rd_offset_kind) :: offset !<
464
465	#if defined( __parallel )
466	INTEGER, DIMENSION(rd_status_size) :: status !<
467	#endif
468
469	LOGICAL, INTENT(IN), OPTIONAL :: open_for_global_io_only !<
470	LOGICAL :: set_filetype !<
471
472	#if ! defined( __parallel )
473	TYPE(C_PTR) :: buf_ptr !<
474	#endif
475
476
477	offset = 0
478	io_on_limited_cores_per_node = .FALSE.
479
480	rd_flag = ( TRIM( action ) == 'READ' .OR. TRIM( action ) == 'read' )
481	wr_flag = ( TRIM( action ) == 'WRITE' .OR. TRIM( action ) == 'write' )
482
483	IF ( .NOT. ( rd_flag .OR. wr_flag ) ) THEN
484	message_string = 'illegal action "' // TRIM( action ) // '" for opening restart files'
485	CALL message( 'restart_data_mpi_io_mod', 'PA0720', 1, 2, 0, 6, 0 )
486	ENDIF
487	!
488	!-- Store name of I/O file to communicate it internally within this module.
489	io_file_name = file_name
490	!
491	!-- Setup for IO on a limited number of PEs per node (using shared memory MPI)
492	IF ( rd_flag ) THEN
493	set_filetype = .TRUE.
494	IF ( TRIM( restart_data_format_input ) == 'mpi_shared_memory' ) THEN
495	io_on_limited_cores_per_node = .TRUE.
496	ENDIF
497	ENDIF
498
499	IF ( TRIM( restart_data_format_output ) == 'mpi_shared_memory' .AND. wr_flag ) THEN
500	io_on_limited_cores_per_node = .TRUE.
501	ENDIF
502	!
503	!-- Shared memory MPI is not used for reading of global data
504	IF ( PRESENT( open_for_global_io_only ) .AND. rd_flag ) THEN
505	IF ( open_for_global_io_only ) THEN
506	io_on_limited_cores_per_node = .FALSE.
507	set_filetype = .FALSE.
508	ENDIF
509	ENDIF
510
511	!
512	!-- Determine, if prerun data shall be read and mapped cyclically to the mainrun arrays.
513	!-- In cyclic fill mode only a subset of the PEs will read.
514	cyclic_fill_mode = .FALSE.
515	pe_active_for_read = .TRUE.
516
517	IF ( rd_flag .AND. .NOT. PRESENT( open_for_global_io_only ) .AND. &
518	nx_on_file < nx .AND. ny_on_file < ny ) THEN
519	cyclic_fill_mode = .TRUE.
520	CALL setup_cyclic_fill
521	!
522	!-- Shared memory IO on limited cores is not allowed for cyclic fill mode
523	CALL sm_io%sm_init_comm( .FALSE. ) !
524	ELSE
525	CALL sm_io%sm_init_comm( io_on_limited_cores_per_node )
526	ENDIF
527
528	!
529	!-- TODO: add a more detailed meaningful comment about what is happening here
530	!-- activate model grid
531	IF( cyclic_fill_mode .AND. .NOT. pe_active_for_read ) THEN
532	CALL mainrun_grid%activate_grid_from_this_class()
533	RETURN
534	ENDIF
535
536
537	!
538	!-- Set communicator to be used. If all cores are doing I/O, comm2d is used as usual.
539	IF( sm_io%is_sm_active() ) THEN
540	comm_io = sm_io%comm_io
541	ELSEIF ( cyclic_fill_mode ) THEN
542	comm_io = comm_cyclic_fill
543	ELSE
544	comm_io = comm2d
545	ENDIF
546
547	!
548	!-- Create subarrays and file types
549	filetypes_created = .FALSE.
550
551	!
552	!-- In case of read it is not known yet if data include total domain. Filetypes will be created
553	!-- further below.
554	IF ( wr_flag ) THEN
555	CALL rd_mpi_io_create_filetypes
556	filetypes_created = .TRUE.
557	ENDIF
558
559	!
560	!-- Open file for MPI-IO
561	#if defined( __parallel )
562	IF ( sm_io%iam_io_pe ) THEN
563
564	IF ( rd_flag ) THEN
565
566	IF ( debug_output ) THEN
567	WRITE( debug_string, * ) 'open joint restart file "' // TRIM( io_file_name ) // &
568	'" for read with MPI-IO'
569	CALL debug_message( debug_string, 'start' )
570	ENDIF
571
572	CALL MPI_FILE_OPEN( comm_io, TRIM( io_file_name ), MPI_MODE_RDONLY, MPI_INFO_NULL, fh, &
573	ierr )
574
575	IF ( ierr /= 0 ) THEN
576	message_string = 'error opening restart file "' // TRIM( io_file_name ) // &
577	'" for reading with MPI-IO'
578	CALL message( 'rrd_mpi_io_open', 'PA0727', 3, 2, 0, 6, 0 )
579	ENDIF
580
581	IF ( debug_output ) THEN
582	WRITE( debug_string, * ) 'open joint restart file "' // TRIM( io_file_name ) // &
583	'" for read with MPI-IO'
584	CALL debug_message( debug_string, 'end' )
585	ENDIF
586
587	ELSEIF ( wr_flag ) THEN
588
589	IF ( debug_output ) THEN
590	WRITE( debug_string, * ) 'open joint restart file "' // TRIM( io_file_name ) // &
591	'" for write with MPI-IO'
592	CALL debug_message( debug_string, 'start' )
593	ENDIF
594
595	CALL MPI_FILE_OPEN( comm_io, TRIM( io_file_name ), MPI_MODE_CREATE+MPI_MODE_WRONLY, &
596	MPI_INFO_NULL, fh, ierr )
597
598	IF ( ierr /= 0 ) THEN
599	message_string = 'error opening restart file "' // TRIM( io_file_name ) // &
600	'" for writing with MPI-IO'
601	CALL message( 'rrd_mpi_io_open', 'PA0728', 3, 2, 0, 6, 0 )
602	ENDIF
603
604	IF ( debug_output ) THEN
605	WRITE( debug_string, * ) 'open joint restart file "' // TRIM( io_file_name ) // &
606	'" for write with MPI-IO'
607	CALL debug_message( debug_string, 'end' )
608	ENDIF
609
610	ENDIF
611
612	ENDIF
613	#else
614	IF ( rd_flag ) THEN
615
616	IF ( debug_output ) THEN
617	WRITE( debug_string, * ) 'open restart file "' // TRIM( io_file_name ) // &
618	'" for read in serial mode (posix)'
619	CALL debug_message( debug_string, 'start' )
620	ENDIF
621
622	fh = posix_open( TRIM( io_file_name ), .TRUE. )
623
624	IF ( debug_output ) THEN
625	WRITE( debug_string, * ) 'open restart file "' // TRIM( io_file_name ) // &
626	'" for read in serial mode (posix)'
627	CALL debug_message( debug_string, 'end' )
628	ENDIF
629
630	ELSEIF ( wr_flag ) THEN
631
632	IF ( debug_output ) THEN
633	WRITE( debug_string, * ) 'open restart file "' // TRIM( io_file_name ) // &
634	'" for write in serial mode (posix)'
635	CALL debug_message( debug_string, 'start' )
636	ENDIF
637
638	fh = posix_open( TRIM( io_file_name ), .FALSE. )
639
640	IF ( debug_output ) THEN
641	WRITE( debug_string, * ) 'open restart file "' // TRIM( io_file_name ) // &
642	'" for write in serial mode (posix)'
643	CALL debug_message( debug_string, 'end' )
644	ENDIF
645
646	ENDIF
647
648	IF ( fh < 0 ) THEN
649	message_string = 'error opening restart file for posix I/O'
650	CALL message( 'restart_data_mpi_io_mod', 'PA0721', 1, 2, 0, 6, 0 )
651	ENDIF
652	#endif
653
654	array_position = 65536 !> Start offset for writing 2-D and 3.D arrays at 64 k
655	header_position = 0
656
657	header_int_index = 1
658	header_char_index = 1
659	header_real_index = 1
660	header_array_index = 1
661
662	int_names = ' '
663	int_values = 0
664	text_lines = ' '
665	real_names = ' '
666	real_values = 0.0
667	array_names = ' '
668	array_offset = 0
669
670	int_names(1) = 'nx'
671	int_values(1) = nx
672	int_names(2) = 'ny'
673	int_values(2) = ny
674	int_names(3) = 'nz'
675	int_values(3) = nz
676	header_int_index = header_int_index+3
677
678	DO i = 1, max_nr_arrays
679	array_offset(i) = 0
680	array_names(i) = ' '
681	ENDDO
682
683	gh_size = STORAGE_SIZE( tgh ) / 8
684
685	IF ( rd_flag ) THEN
686	IF ( sm_io%iam_io_pe ) THEN
687	!
688	!-- File is open for reading
689	#if defined( __parallel )
690	!-- Set the default view
691	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
692	!
693	!-- Read the file header size
694	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
695	CALL MPI_FILE_READ( fh, tgh, gh_size, MPI_BYTE, status, ierr )
696	#else
697	CALL posix_lseek( fh, header_position )
698	buf_ptr = C_LOC( tgh )
699	CALL posix_read( fh, buf_ptr, gh_size )
700	#endif
701	ENDIF
702	#if defined( __parallel )
703	IF ( sm_io%is_sm_active() ) THEN
704	CALL MPI_BCAST( tgh, gh_size, MPI_BYTE, 0, sm_io%comm_shared, ierr )
705	ENDIF
706	#endif
707	header_position = header_position + gh_size
708
709	include_total_domain_boundaries = ( tgh%i_outer_bound == 1 )
710
711	!
712	!-- File types depend on boundaries of the total domain being included in data. This has been
713	!-- checked with the previous statement.
714	IF ( set_filetype ) THEN
715	CALL rd_mpi_io_create_filetypes
716	filetypes_created = .TRUE.
717	ENDIF
718
719	IF ( sm_io%iam_io_pe ) THEN
720	#if defined( __parallel )
721	!
722	!-- Read INTEGER values
723	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
724	CALL MPI_FILE_READ( fh, int_names, SIZE( int_names ) * 32, MPI_CHAR, status, ierr )
725	header_position = header_position + SIZE( int_names ) * 32
726
727	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
728	CALL MPI_FILE_READ (fh, int_values, SIZE( int_values ), MPI_INT, status, ierr )
729	header_position = header_position + SIZE( int_values ) * iwp
730	!
731	!-- Character entries
732	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
733	CALL MPI_FILE_READ( fh, text_lines, SIZE( text_lines ) * 128, MPI_CHAR, status, ierr )
734	header_position = header_position + SIZE ( text_lines ) * 128
735	!
736	!-- REAL values
737	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
738	CALL MPI_FILE_READ( fh, real_names, SIZE( real_names ) * 32, MPI_CHAR, status, ierr )
739	header_position = header_position + SIZE( real_names ) * 32
740
741	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
742	CALL MPI_FILE_READ( fh, real_values, SIZE( real_values ), MPI_REAL, status, ierr )
743	header_position = header_position + SIZE( real_values ) * wp
744	!
745	!-- 2d- and 3d-array headers
746	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
747	CALL MPI_FILE_READ( fh, array_names, SIZE( array_names ) * 32, MPI_CHAR, status, ierr )
748	header_position = header_position + SIZE( array_names ) * 32
749
750	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
751	CALL MPI_FILE_READ( fh, array_offset, SIZE( array_offset ) * MPI_OFFSET_KIND, MPI_BYTE, &
752	status,ierr ) ! there is no I*8 datatype in Fortran
753	header_position = header_position + SIZE( array_offset ) * rd_offset_kind
754	#else
755	CALL posix_lseek( fh, header_position )
756	CALL posix_read( fh, int_names )
757	header_position = header_position + SIZE( int_names ) * 32
758
759	CALL posix_lseek( fh, header_position )
760	CALL posix_read( fh, int_values, SIZE( int_values ) )
761	header_position = header_position + SIZE( int_values ) * iwp
762	!
763	!-- Character entries
764	CALL posix_lseek( fh, header_position )
765	CALL posix_read( fh, text_lines )
766	header_position = header_position + SIZE( text_lines ) * 128
767	!
768	!-- REAL values
769	CALL posix_lseek( fh, header_position )
770	CALL posix_read( fh, real_names )
771	header_position = header_position + SIZE( real_names ) * 32
772
773	CALL posix_lseek( fh, header_position )
774	CALL posix_read( fh, real_values, SIZE( real_values ) )
775	header_position = header_position + SIZE( real_values ) * wp
776	!
777	!-- 2d- and 3d-array headers
778	CALL posix_lseek( fh, header_position )
779	CALL posix_read( fh, array_names )
780	header_position = header_position + SIZE( array_names ) * 32
781
782	CALL posix_lseek( fh, header_position )
783	CALL posix_read( fh, array_offset, SIZE( array_offset ) ) ! there is no I*8 datatype in Fortran
784	header_position = header_position + SIZE( array_offset ) * rd_offset_kind
785	#endif
786	IF ( debug_output ) CALL rd_mpi_io_print_header
787
788	ENDIF
789
790	#if defined( __parallel )
791	!
792	!-- Broadcast header to all remaining cores that are not involved in I/O
793	IF ( sm_io%is_sm_active() ) THEN
794	!
795	!-- Not sure, that it is possible to broadcast CHARACTER array in one MPI_Bcast call
796	DO i = 1, SIZE( int_names )
797	CALL MPI_BCAST( int_names(i), 32, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
798	ENDDO
799	CALL MPI_BCAST( int_values, SIZE( int_values ), MPI_INTEGER, 0, sm_io%comm_shared, ierr )
800
801	DO i = 1, SIZE( text_lines )
802	CALL MPI_BCAST( text_lines(i), 128, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
803	ENDDO
804
805	DO i = 1, SIZE( real_names )
806	CALL MPI_BCAST( real_names(i), 32, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
807	ENDDO
808	CALL MPI_BCAST( real_values, SIZE( real_values ), MPI_REAL, 0, sm_io%comm_shared, ierr )
809
810	DO i = 1, SIZE( array_names )
811	CALL MPI_BCAST( array_names(i), 32, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
812	ENDDO
813	CALL MPI_BCAST( array_offset, SIZE( array_offset )*8, MPI_BYTE, 0, sm_io%comm_shared, &
814	ierr ) ! there is no I8 datatype in Fortran (array_offset is I8!)
815
816	CALL MPI_BCAST( header_position, rd_offset_kind, MPI_BYTE, 0, sm_io%comm_shared, ierr )
817
818	ENDIF
819	#endif
820
821
822	ENDIF
823
824	!
825	!-- TODO: describe in more detail what is happening here
826	!-- activate model grid
827	IF ( cyclic_fill_mode ) CALL mainrun_grid%activate_grid_from_this_class()
828
829	CONTAINS
830
831	SUBROUTINE setup_cyclic_fill
832
833	IMPLICIT NONE
834
835	INTEGER :: color !< used to set the IO PEs for MPI_COMM_SPLIT
836	INTEGER(iwp) :: i !<
837	INTEGER(iwp) :: j !<
838	#if defined( __parallel )
839	INTEGER :: ierr !<
840	INTEGER(KIND=MPI_ADDRESS_KIND) :: winsize !< size of RMA window
841	#else
842	INTEGER(idp) :: winsize
843	#endif
844
845	!
846	!-- Save grid information of the mainrun, i.e. grid variables like nxl, nxr, nys, nyn and other
847	!-- values are stored within the mainrun_grid structure
848	CALL mainrun_grid%save_grid_into_this_class()
849
850	ALLOCATE( remote_pe(0:nx_on_file,0:ny_on_file) )
851	ALLOCATE( remote_pe_s(0:nx_on_file,0:ny_on_file) )
852	ALLOCATE( rma_offset(0:nx_on_file,0:ny_on_file) )
853	ALLOCATE( rma_offset_s(0:nx_on_file,0:ny_on_file) )
854
855	remote_pe_s = 0
856	rma_offset_s = 0
857	!
858	!-- Determine, if gridpoints of the prerun are located on this PE.
859	!-- Set the (cyclic) prerun grid.
860	nxr = MIN( nxr, nx_on_file )
861	IF ( nxl > nx_on_file ) THEN
862	nxl = -99
863	nxr = -99
864	nnx = 0
865	ELSE
866	nnx =nxr-nxl+1
867	ENDIF
868
869	nyn = MIN( nyn, ny_on_file )
870	IF ( nys > ny_on_file ) THEN
871	nys = -99
872	nyn = -99
873	nny = 0
874	ELSE
875	nny = nyn-nys+1
876	ENDIF
877
878	nx = nx_on_file
879	ny = ny_on_file
880	!
881	!-- Determine, if this PE is doing IO
882	IF ( nnx > 0 .AND. nny > 0 ) THEN
883	color = 1
884	pe_active_for_read = .TRUE.
885	remote_pe_s(nxl:nxr,nys:nyn) = myid ! myid from comm2d
886	DO j = nys, nyn
887	DO i = nxl, nxr
888	rma_offset_s(i,j) = ( j-nys ) + ( i-nxl ) * nny
889	ENDDO
890	ENDDO
891	ELSE
892	#if defined( __parallel )
893	color = MPI_UNDEFINED
894	#endif
895	pe_active_for_read = .FALSE.
896	ENDIF
897
898	#if defined( __parallel )
899	CALL MPI_ALLREDUCE( remote_pe_s, remote_pe, SIZE(remote_pe_s), MPI_INTEGER, MPI_SUM, &
900	comm2d, ierr )
901	CALL MPI_ALLREDUCE( rma_offset_s, rma_offset, SIZE(rma_offset_s), MPI_INTEGER, MPI_SUM, &
902	comm2d, ierr )
903	CALL MPI_COMM_SPLIT( comm2d, color, 0, comm_cyclic_fill, ierr )
904
905	IF ( pe_active_for_read ) THEN
906	CALL MPI_COMM_SIZE( comm_cyclic_fill, numprocs, ierr )
907	CALL MPI_COMM_RANK( comm_cyclic_fill, myid, ierr )
908	ENDIF
909	#else
910	remote_pe = remote_pe_s
911	rma_offset = rma_offset_s
912	myid = 0
913	numprocs = 1
914	#endif
915	!
916	!-- Allocate 2d buffers as RMA window, accessible on all PEs
917	IF ( pe_active_for_read ) THEN
918	ALLOCATE( rmabuf_2di(nys:nyn,nxl:nxr) )
919	ELSE
920	ALLOCATE( rmabuf_2di(1,1) )
921	ENDIF
922	winsize = SIZE( rmabuf_2di ) * iwp
923
924	#if defined( __parallel )
925	CALL MPI_WIN_CREATE( rmabuf_2di, winsize, iwp, MPI_INFO_NULL, comm2d, rmawin_2di, ierr )
926	CALL MPI_WIN_FENCE( 0, rmawin_2di, ierr )
927	#endif
928
929	IF ( pe_active_for_read ) THEN
930	ALLOCATE( rmabuf_2d(nys:nyn,nxl:nxr) )
931	ELSE
932	ALLOCATE( rmabuf_2d(1,1) )
933	ENDIF
934	winsize = SIZE( rmabuf_2d ) * wp
935
936	#if defined( __parallel )
937	CALL MPI_WIN_CREATE( rmabuf_2d, winsize, wp, MPI_INFO_NULL, comm2d, rmawin_2d, ierr )
938	CALL MPI_WIN_FENCE( 0, rmawin_2d, ierr )
939	#endif
940
941	!
942	!-- Allocate 3d buffer as RMA window, accessable on all PEs
943	IF ( pe_active_for_read ) THEN
944	ALLOCATE( rmabuf_3d(nzb:nzt+1,nys:nyn,nxl:nxr) )
945	ELSE
946	ALLOCATE( rmabuf_3d(1,1,1) )
947	ENDIF
948	winsize = SIZE( rmabuf_3d ) * wp
949
950	#if defined( __parallel )
951	CALL MPI_WIN_CREATE( rmabuf_3d, winsize, wp, MPI_INFO_NULL, comm2d, rmawin_3d, ierr )
952	CALL MPI_WIN_FENCE( 0, rmawin_3d, ierr )
953	#endif
954
955	!
956	!-- Save grid of the prerun, i.e. grid variables like nxl, nxr, nys, nyn and other values
957	!-- are stored within the prerun_grid structure.
958	!-- The prerun grid can later be activated by calling prerun_grid%activate_grid_from_this_class()
959	CALL prerun_grid%save_grid_into_this_class()
960	prerun_grid%comm2d = comm_cyclic_fill
961
962	DEALLOCATE( remote_pe_s, rma_offset_s )
963
964	END SUBROUTINE setup_cyclic_fill
965
966	END SUBROUTINE rd_mpi_io_open
967
968
969	!--------------------------------------------------------------------------------------------------!
970	! Description:
971	! ------------
972	!> Check, if array exists in restart file
973	!--------------------------------------------------------------------------------------------------!
974	SUBROUTINE rd_mpi_io_check_array( name, found )
975
976	IMPLICIT NONE
977
978	CHARACTER(LEN=*), INTENT(IN) :: name !<
979
980	INTEGER(iwp) :: i !<
981
982	LOGICAl :: found !<
983
984
985	DO i = 1, tgh%nr_arrays
986	IF ( TRIM( array_names(i) ) == TRIM( name ) ) THEN
987	array_position = array_offset(i)
988	found = .TRUE.
989	RETURN
990	ENDIF
991	ENDDO
992
993	found = .FALSE.
994
995	END SUBROUTINE rd_mpi_io_check_array
996
997
998
999	!--------------------------------------------------------------------------------------------------!
1000	! Description:
1001	! ------------
1002	!> Read INTEGER with MPI-IO
1003	!--------------------------------------------------------------------------------------------------!
1004	SUBROUTINE rrd_mpi_io_int( name, value )
1005
1006	IMPLICIT NONE
1007
1008	CHARACTER(LEN=*), INTENT(IN) :: name !<
1009
1010	INTEGER(iwp) :: i !<
1011	INTEGER(KIND=iwp), INTENT(OUT) :: value !<
1012
1013	LOGICAL :: found !<
1014
1015
1016	found = .FALSE.
1017	value = 0
1018
1019	DO i = 1, tgh%nr_int
1020	IF ( TRIM(int_names(i)) == TRIM( name ) ) THEN
1021	value = int_values(i)
1022	found = .TRUE.
1023	EXIT
1024	ENDIF
1025	ENDDO
1026
1027	IF ( .NOT. found ) THEN
1028	message_string = 'INTEGER variable "' // TRIM( name ) // '" not found in restart file'
1029	CALL message( 'rrd_mpi_io_int', 'PA0722', 3, 2, 0, 6, 0 )
1030	ENDIF
1031
1032	END SUBROUTINE rrd_mpi_io_int
1033
1034
1035
1036	!--------------------------------------------------------------------------------------------------!
1037	! Description:
1038	! ------------
1039	!> Read REAL with MPI-IO
1040	!--------------------------------------------------------------------------------------------------!
1041	SUBROUTINE rrd_mpi_io_real( name, value )
1042
1043	IMPLICIT NONE
1044
1045	CHARACTER(LEN=*), INTENT(IN) :: name !<
1046
1047	INTEGER(iwp) :: i !<
1048
1049	LOGICAL :: found !<
1050
1051	REAL(KIND=wp), INTENT(OUT) :: value !<
1052
1053
1054	found = .FALSE.
1055	value = 0.0
1056
1057	DO i = 1, tgh%nr_real
1058	IF ( TRIM(real_names(i)) == TRIM( name ) ) THEN
1059	value = real_values(i)
1060	found = .TRUE.
1061	EXIT
1062	ENDIF
1063	ENDDO
1064
1065	IF ( .NOT. found ) THEN
1066	message_string = 'REAL variable "' // TRIM( name ) // '" not found in restart file'
1067	CALL message( 'rrd_mpi_io_int', 'PA0722', 3, 2, 0, 6, 0 )
1068	ENDIF
1069
1070	END SUBROUTINE rrd_mpi_io_real
1071
1072
1073
1074	!--------------------------------------------------------------------------------------------------!
1075	! Description:
1076	! ------------
1077	!> Read 2d-real array with MPI-IO
1078	!--------------------------------------------------------------------------------------------------!
1079	SUBROUTINE rrd_mpi_io_real_2d( name, data )
1080
1081	IMPLICIT NONE
1082
1083	CHARACTER(LEN=*), INTENT(IN) :: name !<
1084
1085	#if defined( __parallel )
1086	INTEGER, DIMENSION(rd_status_size) :: status !<
1087	#endif
1088	INTEGER(iwp) :: i !<
1089
1090	LOGICAL :: found !<
1091
1092	REAL(wp), INTENT(INOUT), DIMENSION(nysg:nyng,nxlg:nxrg) :: data !<
1093
1094
1095	found = .FALSE.
1096
1097	DO i = 1, tgh%nr_arrays
1098	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
1099	array_position = array_offset(i)
1100	found = .TRUE.
1101	EXIT
1102	ENDIF
1103	ENDDO
1104
1105	IF ( found ) THEN
1106
1107	IF ( cyclic_fill_mode ) THEN
1108
1109	CALL rrd_mpi_io_real_2d_cyclic_fill
1110
1111	ELSE
1112
1113	#if defined( __parallel )
1114	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited # of cores is inactive
1115	IF ( sm_io%iam_io_pe ) THEN
1116	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL, &
1117	ierr )
1118	CALL MPI_FILE_READ_ALL( fh, array_2d, SIZE( array_2d ), MPI_REAL, status, ierr )
1119	ENDIF
1120	CALL sm_io%sm_node_barrier()
1121	#else
1122	CALL posix_lseek( fh, array_position )
1123	CALL posix_read( fh, array_2d, SIZE( array_2d ) )
1124	#endif
1125
1126	IF ( include_total_domain_boundaries ) THEN
1127	DO i = iog%nxl, iog%nxr
1128	data(iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_2d(i,iog%nys:iog%nyn)
1129	ENDDO
1130	IF ( debug_level >= 2) THEN
1131	WRITE(9,*) 'r2f_ob ', TRIM(name),' ', SUM( data(nys:nyn,nxl:nxr) )
1132	ENDIF
1133	ELSE
1134	DO i = nxl, nxr
1135	data(nys:nyn,i) = array_2d(i,nys:nyn)
1136	ENDDO
1137	IF ( debug_level >= 2) THEN
1138	WRITE(9,*) 'r2f ', TRIM( name ),' ', SUM( data(nys:nyn,nxl:nxr) )
1139	ENDIF
1140	ENDIF
1141
1142	ENDIF
1143
1144	CALL exchange_horiz_2d( data )
1145
1146	ELSE
1147	message_string = '2d-REAL array "' // TRIM( name ) // '" not found in restart file'
1148	CALL message( 'rrd_mpi_io_int', 'PA0722', 3, 2, 0, 6, 0 )
1149	ENDIF
1150
1151
1152	CONTAINS
1153
1154	SUBROUTINE rrd_mpi_io_real_2d_cyclic_fill
1155
1156	IMPLICIT NONE
1157
1158	INTEGER(iwp) :: i !<
1159	INTEGER(iwp) :: ie !<
1160	INTEGER(iwp) :: is !<
1161	INTEGER(iwp) :: i_remote !<
1162	INTEGER(iwp) :: j !<
1163	INTEGER(iwp) :: je !<
1164	INTEGER(iwp) :: js !<
1165	INTEGER(iwp) :: j_remote !<
1166	INTEGER(iwp) :: nval !<
1167	INTEGER(iwp) :: rem_pe !<
1168
1169	#if defined( __parallel )
1170	INTEGER(iwp) :: ierr !<
1171	INTEGER(KIND=MPI_ADDRESS_KIND) :: rem_offs !<
1172	#else
1173	INTEGER(idp) :: rem_offs
1174	#endif
1175
1176
1177	!
1178	!-- Reading 2d real array on prerun grid
1179	CALL prerun_grid%activate_grid_from_this_class()
1180
1181	IF ( pe_active_for_read ) THEN
1182	#if defined( __parallel )
1183	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL, &
1184	ierr )
1185	CALL MPI_FILE_READ_ALL( fh, array_2d, SIZE( array_2d ), MPI_REAL, status, ierr )
1186	#endif
1187	DO i = nxl, nxr
1188	rmabuf_2d(nys:nyn,i) = array_2d(i,nys:nyn)
1189	ENDDO
1190	data(nys:nyn,nxl:nxr) = rmabuf_2d ! copy prerund data directly into output array data
1191	ENDIF
1192
1193	CALL mainrun_grid%activate_grid_from_this_class()
1194
1195	#if defined( __parallel )
1196	!
1197	!-- Close RMA window to allow remote access
1198	CALL MPI_WIN_FENCE( 0, rmawin_2d, ierr )
1199	#endif
1200
1201	!
1202	!-- TODO: describe in more detail what is happening in this IF/ELSE clause
1203	IF ( .NOT. pe_active_for_read ) THEN
1204
1205	is = nxl
1206	ie = nxr
1207	js = nys
1208	je = nyn
1209
1210	ELSE
1211	!
1212	!-- Extra get for cyclic data north of prerun data
1213	is = nxl
1214	ie = nxr
1215	js = prerun_grid%nys+1
1216	je = nyn
1217	DO i = is, ie
1218	DO j = js, je
1219	i_remote = MOD(i,nx_on_file+1)
1220	j_remote = MOD(j,ny_on_file+1)
1221	rem_pe = remote_pe(i_remote,j_remote)
1222	rem_offs = rma_offset(i_remote,j_remote)
1223	nval = 1
1224
1225	#if defined( __parallel )
1226	IF ( rem_pe /= myid ) THEN
1227	CALL MPI_GET( data(j,i), nval, MPI_REAL, rem_pe, rem_offs, nval, MPI_REAL, &
1228	rmawin_2d, ierr )
1229	ELSE
1230	data(j,i) = rmabuf_2d(j_remote,i_remote)
1231	ENDIF
1232	#else
1233	data(j,i) = array_2d(i_remote,j_remote)
1234	#endif
1235	ENDDO
1236	ENDDO
1237	!
1238	!-- Prepare setup for stripe right of prerun data
1239	is = prerun_grid%nxr+1
1240	ie = nxr
1241	js = nys
1242	je = nyn
1243
1244	ENDIF
1245
1246	DO i = is, ie
1247	DO j = js, je
1248	i_remote = MOD(i,nx_on_file+1)
1249	j_remote = MOD(j,ny_on_file+1)
1250	rem_pe = remote_pe(i_remote,j_remote)
1251	rem_offs = rma_offset(i_remote,j_remote)
1252	nval = 1
1253
1254	#if defined( __parallel )
1255	IF ( rem_pe /= myid ) THEN
1256	CALL MPI_GET( data(j,i), nval, MPI_REAL, rem_pe, rem_offs, nval, MPI_REAL, &
1257	rmawin_2d, ierr )
1258	ELSE
1259	data(j,i) = rmabuf_2d(j_remote,i_remote)
1260	ENDIF
1261	#else
1262	data(j,i) = array_2d(i_remote,j_remote)
1263	#endif
1264	ENDDO
1265	ENDDO
1266
1267	#if defined( __parallel )
1268	!
1269	!-- Reopen RMA window to allow filling
1270	CALL MPI_WIN_FENCE( 0, rmawin_2d, ierr )
1271	#endif
1272
1273	END SUBROUTINE rrd_mpi_io_real_2d_cyclic_fill
1274
1275	END SUBROUTINE rrd_mpi_io_real_2d
1276
1277
1278
1279	!--------------------------------------------------------------------------------------------------!
1280	! Description:
1281	! ------------
1282	!> Read 2d-INTEGER array with MPI-IO
1283	!--------------------------------------------------------------------------------------------------!
1284	SUBROUTINE rrd_mpi_io_int_2d( name, data )
1285
1286	IMPLICIT NONE
1287
1288	CHARACTER(LEN=*), INTENT(IN) :: name !<
1289
1290	INTEGER(iwp) :: i !<
1291	INTEGER(iwp) :: j !<
1292
1293	#if defined( __parallel )
1294	INTEGER, DIMENSION(rd_status_size) :: status !<
1295	#endif
1296
1297	INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:,:) :: data !<
1298
1299	LOGICAL :: found !<
1300
1301
1302	found = .FALSE.
1303
1304	DO i = 1, tgh%nr_arrays
1305	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
1306	array_position = array_offset(i)
1307	found = .TRUE.
1308	EXIT
1309	ENDIF
1310	ENDDO
1311
1312	IF ( found ) THEN
1313
1314	IF ( ( nxr - nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) ) THEN
1315	!
1316	!-- Output array with Halos.
1317	!-- ATTENTION: INTEGER arrays with ghost boundaries are not implemented yet. This kind of
1318	!-- array would be dimensioned in the caller subroutine like this:
1319	!-- INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg):: data
1320	message_string = '2d-INTEGER array with nbgp "' // TRIM( name ) // '" to be read ' // &
1321	'from restart file is defined with illegal dimensions in the PALM code'
1322	CALL message( 'rrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
1323
1324	ELSEIF ( (nxr - nxl + 1) == SIZE( data, 2 ) ) THEN
1325	!
1326	!-- INTEGER input array without Halos.
1327	!-- This kind of array is dimensioned in the caller subroutine
1328	!-- INTEGER, DIMENSION(nys:nyn,nxl:nxr) :: data
1329	IF ( cyclic_fill_mode ) THEN
1330
1331	CALL rrd_mpi_io_int_2d_cyclic_fill
1332
1333	ELSE
1334
1335	#if defined( __parallel )
1336	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited # of cores is inactive
1337	IF ( sm_io%iam_io_pe ) THEN
1338	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native', &
1339	MPI_INFO_NULL, ierr )
1340	CALL MPI_FILE_READ_ALL( fh, array_2di, SIZE( array_2di ), MPI_INTEGER, status, &
1341	ierr )
1342	ENDIF
1343	CALL sm_io%sm_node_barrier()
1344	#else
1345	CALL posix_lseek( fh, array_position )
1346	CALL posix_read( fh, array_2di, SIZE( array_2di ) )
1347	#endif
1348	DO j = nys, nyn
1349	DO i = nxl, nxr
1350	data(j-nys+1,i-nxl+1) = array_2di(i,j)
1351	ENDDO
1352	ENDDO
1353
1354	ENDIF
1355
1356	ELSE
1357
1358	message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be read from restart ' // &
1359	'file is defined with illegal dimensions in the PALM code'
1360	CALL message( 'rrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
1361
1362	ENDIF
1363
1364	ELSE
1365
1366	message_string = '2d-INTEGER array "' // TRIM( name ) // '" not found in restart file'
1367	CALL message( 'rrd_mpi_io_int_2d', 'PA0722', 3, 2, 0, 6, 0 )
1368
1369	ENDIF
1370
1371
1372	CONTAINS
1373
1374	SUBROUTINE rrd_mpi_io_int_2d_cyclic_fill
1375
1376	IMPLICIT NONE
1377
1378	INTEGER(iwp) :: i !<
1379	INTEGER(iwp) :: ie !<
1380	INTEGER(iwp) :: is !<
1381	INTEGER(iwp) :: i_remote !<
1382	INTEGER(iwp) :: j !<
1383	INTEGER(iwp) :: je !<
1384	INTEGER(iwp) :: js !<
1385	INTEGER(iwp) :: j_remote !<
1386	INTEGER(iwp) :: nval !<
1387	INTEGER(iwp) :: rem_pe !<
1388
1389	#if defined( __parallel )
1390	INTEGER(iwp) :: ierr !<
1391	INTEGER(KIND=MPI_ADDRESS_KIND) :: rem_offs !<
1392	#else
1393	INTEGER(idp) :: rem_offs
1394	#endif
1395
1396
1397	IF ( pe_active_for_read ) THEN
1398	CALL prerun_grid%activate_grid_from_this_class()
1399
1400	#if defined( __parallel )
1401	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native', &
1402	MPI_INFO_NULL, ierr )
1403	CALL MPI_FILE_READ_ALL( fh, array_2di, SIZE( array_2di ), MPI_INTEGER, status, ierr )
1404	#else
1405	CALL posix_lseek( fh, array_position )
1406	CALL posix_read( fh, array_2di, SIZE( array_2di ) )
1407	#endif
1408	DO i = nxl, nxr
1409	rmabuf_2di(nys:nyn,i) = array_2di(i,nys:nyn)
1410	ENDDO
1411	data(1:nny,1:nnx) = rmabuf_2di
1412
1413	CALL mainrun_grid%activate_grid_from_this_class()
1414	ENDIF
1415
1416	#if defined( __parallel )
1417	!
1418	!-- Close RMA window to allow remote access
1419	CALL MPI_WIN_FENCE( 0, rmawin_2di, ierr )
1420	#endif
1421
1422	is = nxl
1423	ie = nxr
1424	js = nys
1425	je = nyn
1426
1427	DO i = is, ie
1428	DO j = js, je
1429	i_remote = MOD(i,nx_on_file+1)
1430	j_remote = MOD(j,ny_on_file+1)
1431	rem_pe = remote_pe(i_remote,j_remote)
1432	rem_offs = rma_offset(i_remote,j_remote)
1433	nval = 1
1434	#if defined( __parallel )
1435	IF ( rem_pe /= myid ) THEN
1436	CALL MPI_GET( data(j-nys+1,i-nxl+1), nval, MPI_INTEGER, rem_pe, rem_offs, nval, &
1437	MPI_INTEGER, rmawin_2di, ierr)
1438	ELSE
1439	data(j-nys+1,i-nxl+1) = rmabuf_2di(j_remote,i_remote)
1440	ENDIF
1441	#else
1442	data(j-nys+1,i-nxl+1) = array_2di(i_remote,j_remote)
1443	#endif
1444	ENDDO
1445	ENDDO
1446
1447	#if defined( __parallel )
1448	!
1449	!-- Reopen RMA window to allow filling
1450	CALL MPI_WIN_FENCE( 0, rmawin_2di, ierr )
1451	#endif
1452
1453	END SUBROUTINE rrd_mpi_io_int_2d_cyclic_fill
1454
1455	END SUBROUTINE rrd_mpi_io_int_2d
1456
1457
1458
1459	!--------------------------------------------------------------------------------------------------!
1460	! Description:
1461	! ------------
1462	!> Read 3d-INTEGER*4 array with MPI-IO
1463	!--------------------------------------------------------------------------------------------------!
1464	SUBROUTINE rrd_mpi_io_int4_3d( name, data )
1465
1466	IMPLICIT NONE
1467
1468	CHARACTER(LEN=*), INTENT(IN) :: name !<
1469
1470	INTEGER(iwp) :: i !<
1471
1472	#if defined( __parallel )
1473	INTEGER, DIMENSION(rd_status_size) :: status !<
1474	#endif
1475
1476	LOGICAL :: found !<
1477
1478	INTEGER(isp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
1479
1480
1481	found = .FALSE.
1482	data = -1.0
1483
1484	DO i = 1, tgh%nr_arrays
1485	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
1486	array_position = array_offset(i)
1487	found = .TRUE.
1488	EXIT
1489	ENDIF
1490	ENDDO
1491
1492	IF ( found ) THEN
1493
1494	#if defined( __parallel )
1495	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited # of cores is inactive
1496	IF( sm_io%iam_io_pe ) THEN
1497	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_3di4, 'native', &
1498	MPI_INFO_NULL, ierr )
1499	CALL MPI_FILE_READ_ALL( fh, array_3di4, SIZE( array_3di4 ), MPI_INTEGER, status, ierr )
1500	ENDIF
1501	CALL sm_io%sm_node_barrier()
1502	#else
1503	CALL posix_lseek( fh, array_position )
1504	CALL posix_read(fh, array_3di4, SIZE( array_3di4 ) )
1505	#endif
1506	IF ( include_total_domain_boundaries ) THEN
1507	DO i = iog%nxl, iog%nxr
1508	data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3di4(:,i,iog%nys:iog%nyn)
1509	ENDDO
1510	ELSE
1511	DO i = nxl, nxr
1512	data(:,nys:nyn,i) = array_3di4(:,i,nys:nyn)
1513	ENDDO
1514	ENDIF
1515
1516	ELSE
1517
1518	message_string = '3d-INTEGER*4 array "' // TRIM( name ) // '" not found in restart file'
1519	CALL message( 'rrd_mpi_io_int4_3d', 'PA0722', 3, 2, 0, 6, 0 )
1520
1521	ENDIF
1522
1523
1524	END SUBROUTINE rrd_mpi_io_int4_3d
1525
1526
1527
1528	!--------------------------------------------------------------------------------------------------!
1529	! Description:
1530	! ------------
1531	!> Read 3d-INTEGER*8 array with MPI-IO
1532	!--------------------------------------------------------------------------------------------------!
1533	SUBROUTINE rrd_mpi_io_int8_3d( name, data )
1534
1535	IMPLICIT NONE
1536
1537	CHARACTER(LEN=*), INTENT(IN) :: name !<
1538
1539	INTEGER(iwp) :: i !<
1540
1541	#if defined( __parallel )
1542	INTEGER, DIMENSION(rd_status_size) :: status !<
1543	#endif
1544
1545	LOGICAL :: found !<
1546
1547	INTEGER(idp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
1548
1549
1550	found = .FALSE.
1551	data = -1.0
1552
1553	DO i = 1, tgh%nr_arrays
1554	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
1555	array_position = array_offset(i)
1556	found = .TRUE.
1557	EXIT
1558	ENDIF
1559	ENDDO
1560
1561	IF ( found ) THEN
1562
1563	#if defined( __parallel )
1564	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited # of cores is inactive
1565	IF( sm_io%iam_io_pe ) THEN
1566	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER8, ft_3di8, 'native', MPI_INFO_NULL, &
1567	ierr )
1568	CALL MPI_FILE_READ_ALL( fh, array_3di8, SIZE( array_3di8 ), MPI_INTEGER8, status, ierr )
1569	ENDIF
1570	CALL sm_io%sm_node_barrier()
1571	#else
1572	CALL posix_lseek( fh, array_position )
1573	CALL posix_read(fh, array_3di8, SIZE( array_3di8 ) )
1574	#endif
1575	IF ( include_total_domain_boundaries ) THEN
1576	DO i = iog%nxl, iog%nxr
1577	data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3di8(:,i,iog%nys:iog%nyn)
1578	ENDDO
1579	ELSE
1580	DO i = nxl, nxr
1581	data(:,nys:nyn,i) = array_3di8(:,i,nys:nyn)
1582	ENDDO
1583	ENDIF
1584
1585	ELSE
1586
1587	message_string = '3d-INTEGER*8 array "' // TRIM( name ) // '" not found in restart file'
1588	CALL message( 'rrd_mpi_io_int8_3d', 'PA0722', 3, 2, 0, 6, 0 )
1589
1590	ENDIF
1591
1592
1593	END SUBROUTINE rrd_mpi_io_int8_3d
1594
1595
1596
1597	!--------------------------------------------------------------------------------------------------!
1598	! Description:
1599	! ------------
1600	!> Read 2d-REAL array with MPI-IO
1601	!--------------------------------------------------------------------------------------------------!
1602	SUBROUTINE rrd_mpi_io_real_3d( name, data )
1603
1604	IMPLICIT NONE
1605
1606	CHARACTER(LEN=*), INTENT(IN) :: name !<
1607
1608	INTEGER(iwp) :: i !<
1609
1610	#if defined( __parallel )
1611	INTEGER, DIMENSION(rd_status_size) :: status !<
1612	#endif
1613
1614	LOGICAL :: found !<
1615
1616	REAL(wp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
1617
1618
1619	found = .FALSE.
1620	data = -1.0
1621
1622	DO i = 1, tgh%nr_arrays
1623	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
1624	array_position = array_offset(i)
1625	found = .TRUE.
1626	EXIT
1627	ENDIF
1628	ENDDO
1629
1630	IF ( found ) THEN
1631
1632	IF ( cyclic_fill_mode ) THEN
1633
1634	CALL rrd_mpi_io_real_3d_cyclic_fill
1635	!
1636	!-- Cyclic fill mode requires to use the "cyclic" communicator, in order to initialize
1637	!-- grid points at the outer boundaries (ghost layers) of the total domain. These points
1638	!-- are not contained in the prerun data, because the prerun used cyclic boundary conditions.
1639	CALL exchange_horiz( data, nbgp, alternative_communicator = 1 )
1640
1641	ELSE
1642	#if defined( __parallel )
1643	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited # of cores is inactive
1644	IF( sm_io%iam_io_pe ) THEN
1645	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL, &
1646	ierr )
1647	CALL MPI_FILE_READ_ALL( fh, array_3d, SIZE( array_3d ), MPI_REAL, status, ierr )
1648	ENDIF
1649	CALL sm_io%sm_node_barrier()
1650	#else
1651	CALL posix_lseek( fh, array_position )
1652	CALL posix_read(fh, array_3d, SIZE( array_3d ) )
1653	#endif
1654	IF ( include_total_domain_boundaries ) THEN
1655	DO i = iog%nxl, iog%nxr
1656	data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3d(:,i,iog%nys:iog%nyn)
1657	ENDDO
1658	ELSE
1659	DO i = nxl, nxr
1660	data(:,nys:nyn,i) = array_3d(:,i,nys:nyn)
1661	ENDDO
1662	ENDIF
1663
1664	CALL exchange_horiz( data, nbgp )
1665
1666	ENDIF
1667
1668	ELSE
1669
1670	message_string = '3d-REAL array "' // TRIM( name ) // '" not found in restart file'
1671	CALL message( 'rrd_mpi_io_real_3d', 'PA0722', 3, 2, 0, 6, 0 )
1672
1673	ENDIF
1674
1675
1676	CONTAINS
1677
1678	SUBROUTINE rrd_mpi_io_real_3d_cyclic_fill
1679
1680	IMPLICIT NONE
1681
1682	INTEGER(iwp) :: i !<
1683	INTEGER(iwp) :: ie !<
1684	INTEGER(iwp) :: is !<
1685	INTEGER(iwp) :: i_remote !<
1686	INTEGER(iwp) :: j !<
1687	INTEGER(iwp) :: je !<
1688	INTEGER(iwp) :: js !<
1689	INTEGER(iwp) :: j_remote !<
1690	INTEGER(iwp) :: nval !<
1691	INTEGER(iwp) :: rem_pe !<
1692
1693	#if defined( __parallel )
1694	INTEGER(iwp) :: ierr !<
1695	INTEGER(KIND=MPI_ADDRESS_KIND) :: rem_offs !<
1696	#else
1697	INTEGER(idp) :: rem_offs
1698	#endif
1699
1700
1701	CALL prerun_grid%activate_grid_from_this_class()
1702
1703	IF ( pe_active_for_read ) THEN
1704	#if defined( __parallel )
1705	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL, &
1706	ierr )
1707	CALL MPI_FILE_READ_ALL( fh, array_3d, SIZE( array_3d ), MPI_REAL, status, ierr )
1708	#else
1709	CALL posix_lseek( fh, array_position )
1710	CALL posix_read( fh, array_3d, SIZE( array_3d ) )
1711	#endif
1712	DO i = nxl, nxr
1713	rmabuf_3d(:,nys:nyn,i) = array_3d(:,i,nys:nyn)
1714	ENDDO
1715	data(:,nys:nyn,nxl:nxr) = rmabuf_3d
1716	ENDIF
1717	CALL mainrun_grid%activate_grid_from_this_class ()
1718
1719	#if defined( __parallel )
1720	!
1721	!-- Close RMA window to allow remote access
1722	CALL MPI_WIN_FENCE( 0, rmawin_3d, ierr )
1723	#endif
1724
1725	is = nxl
1726	ie = nxr
1727	js = nys
1728	je = nyn
1729
1730	DO i = is, ie
1731	DO j = js, je
1732	i_remote = MOD( i, nx_on_file+1 )
1733	j_remote = MOD( j, ny_on_file+1 )
1734	rem_pe = remote_pe(i_remote,j_remote)
1735	rem_offs = rma_offset(i_remote,j_remote) * ( nzt-nzb+2 )
1736	nval = nzt-nzb+2
1737
1738	#if defined( __parallel )
1739	IF ( rem_pe /= myid ) THEN
1740	CALL MPI_GET( data(nzb,j,i), nval, MPI_REAL, rem_pe, rem_offs, nval, MPI_REAL, &
1741	rmawin_3d, ierr)
1742	ELSE
1743	data(:,j,i) = rmabuf_3d(:,j_remote,i_remote)
1744	ENDIF
1745	#else
1746	data(:,j,i) = array_3d(:,i_remote,j_remote)
1747	#endif
1748	ENDDO
1749	ENDDO
1750
1751	#if defined( __parallel )
1752	!
1753	!-- Reopen RMA window to allow filling
1754	CALL MPI_WIN_FENCE( 0, rmawin_3d, ierr )
1755	#endif
1756
1757	END SUBROUTINE rrd_mpi_io_real_3d_cyclic_fill
1758
1759	END SUBROUTINE rrd_mpi_io_real_3d
1760
1761
1762
1763	!--------------------------------------------------------------------------------------------------!
1764	! Description:
1765	! ------------
1766	!> Read 3d-REAL soil array with MPI-IO
1767	!> nzb_soil, nzt_soil are located in the module land_surface_model_mod. Since Fortran does not allow
1768	!> cross referencing of module variables, it is required to pass these variables as arguments.
1769	!--------------------------------------------------------------------------------------------------!
1770	SUBROUTINE rrd_mpi_io_real_3d_soil( name, data, nzb_soil, nzt_soil )
1771
1772	IMPLICIT NONE
1773
1774	CHARACTER(LEN=*), INTENT(IN) :: name !<
1775
1776	INTEGER(iwp) :: i !<
1777	INTEGER, INTENT(IN) :: nzb_soil !<
1778	INTEGER, INTENT(IN) :: nzt_soil !<
1779
1780	#if defined( __parallel )
1781	INTEGER, DIMENSION(rd_status_size) :: status !<
1782	INTEGER(iwp) :: ierr !<
1783	#endif
1784
1785	LOGICAL :: found !<
1786
1787	REAL(wp), INTENT(INOUT), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) :: data !<
1788
1789
1790	!
1791	!-- Prerun data is not allowed to contain soil information so far
1792	IF ( cyclic_fill_mode ) THEN
1793	message_string = 'prerun data is not allowed to contain soil information'
1794	CALL message( 'rrd_mpi_io_real_3d_soil', 'PA0729', 3, 2, -1, 6, 0 )
1795	ENDIF
1796
1797	found = .FALSE.
1798
1799	DO i = 1, tgh%nr_arrays
1800	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
1801	array_position = array_offset(i)
1802	found = .TRUE.
1803	EXIT
1804	ENDIF
1805	ENDDO
1806
1807	IF ( found ) THEN
1808	CALL rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
1809	#if defined( __parallel )
1810	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
1811	IF ( sm_io%iam_io_pe ) THEN
1812	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native', &
1813	MPI_INFO_NULL, ierr )
1814	CALL MPI_FILE_READ_ALL( fh, array_3d_soil, SIZE( array_3d_soil ), MPI_REAL, status, ierr )
1815	CALL MPI_TYPE_FREE( ft_3dsoil, ierr )
1816	ENDIF
1817	CALL sm_io%sm_node_barrier()
1818	#else
1819	CALL posix_lseek( fh, array_position )
1820	CALL posix_read( fh, array_3d_soil, SIZE( array_3d_soil ) )
1821	#endif
1822	IF ( include_total_domain_boundaries ) THEN
1823	DO i = iog%nxl, iog%nxr
1824	data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3d_soil(:,i,iog%nys:iog%nyn)
1825	ENDDO
1826	ELSE
1827	DO i = nxl, nxr
1828	data(:,nys:nyn,i) = array_3d_soil(:,i,nys:nyn)
1829	ENDDO
1830	ENDIF
1831
1832	#if defined( __parallel )
1833	IF ( sm_io%is_sm_active() ) THEN
1834	CALL MPI_WIN_FREE( win_3ds, ierr )
1835	ELSE
1836	DEALLOCATE( array_3d_soil )
1837	ENDIF
1838	#else
1839	DEALLOCATE( array_3d_soil )
1840	#endif
1841
1842	ELSE
1843
1844	message_string = '3d-REAL soil array "' // TRIM( name ) // '" not found in restart file'
1845	CALL message( 'rrd_mpi_io_real_3d_soil', 'PA0722', 3, 2, 0, 6, 0 )
1846
1847	ENDIF
1848
1849	END SUBROUTINE rrd_mpi_io_real_3d_soil
1850
1851
1852
1853	!--------------------------------------------------------------------------------------------------!
1854	! Description:
1855	! ------------
1856	!> Read CHARACTER with MPI-IO
1857	!--------------------------------------------------------------------------------------------------!
1858	SUBROUTINE rrd_mpi_io_char( name, text )
1859
1860	IMPLICIT NONE
1861
1862	CHARACTER(LEN=*), INTENT(IN) :: name !<
1863	CHARACTER(LEN=*), INTENT(OUT) :: text !<
1864	CHARACTER(LEN=128) :: line !<
1865
1866	INTEGER(iwp) :: i !<
1867
1868	LOGICAL :: found !<
1869
1870
1871	found = .FALSE.
1872	text = ' '
1873
1874	DO i = 1, tgh%nr_char
1875	line = text_lines(i)
1876	IF ( TRIM( line(1:32) ) == TRIM( name ) ) THEN
1877	text = line(33:)
1878	found = .TRUE.
1879	EXIT
1880	ENDIF
1881	ENDDO
1882
1883	IF ( .NOT. found ) THEN
1884	message_string = 'CHARACTER variable "' // TRIM( name ) // '" not found in restart file'
1885	CALL message( 'rrd_mpi_io_char', 'PA0722', 3, 2, 0, 6, 0 )
1886	ENDIF
1887
1888	END SUBROUTINE rrd_mpi_io_char
1889
1890
1891
1892	!--------------------------------------------------------------------------------------------------!
1893	! Description:
1894	! ------------
1895	!> Read LOGICAL with MPI-IO
1896	!--------------------------------------------------------------------------------------------------!
1897	SUBROUTINE rrd_mpi_io_logical( name, value )
1898
1899	IMPLICIT NONE
1900
1901	CHARACTER(LEN=*), INTENT(IN) :: name !<
1902
1903	INTEGER(iwp) :: logical_as_integer !<
1904
1905	LOGICAL, INTENT(OUT) :: value !<
1906
1907
1908	CALL rrd_mpi_io_int( name, logical_as_integer )
1909	value = ( logical_as_integer == 1 )
1910
1911	END SUBROUTINE rrd_mpi_io_logical
1912
1913
1914
1915	!--------------------------------------------------------------------------------------------------!
1916	! Description:
1917	! ------------
1918	!> Read 3D-LOGICAL with MPI-IO
1919	!--------------------------------------------------------------------------------------------------!
1920	SUBROUTINE rrd_mpi_io_logical_3d( name, data )
1921
1922	IMPLICIT NONE
1923
1924	CHARACTER(LEN=*), INTENT(IN) :: name !<
1925
1926	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: logical_as_integer !<
1927
1928	LOGICAL, INTENT(OUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
1929
1930
1931	CALL rrd_mpi_io_int4_3d( name, logical_as_integer )
1932	data(:,:,:) = ( logical_as_integer(:,:,:) == 1 )
1933
1934	END SUBROUTINE rrd_mpi_io_logical_3d
1935
1936
1937
1938	!--------------------------------------------------------------------------------------------------!
1939	! Description:
1940	! ------------
1941	!> Write INTEGER with MPI-IO
1942	!--------------------------------------------------------------------------------------------------!
1943	SUBROUTINE wrd_mpi_io_int( name, value )
1944
1945	IMPLICIT NONE
1946
1947	CHARACTER(LEN=*), INTENT(IN) :: name !<
1948
1949	INTEGER(KIND=iwp), INTENT(IN) :: value !<
1950
1951
1952	IF ( header_int_index == max_nr_int ) THEN
1953	STOP '+++ maximum number of INTEGER entries in restart file header exceeded'
1954	ENDIF
1955
1956	int_names(header_int_index) = name
1957	int_values(header_int_index) = value
1958	header_int_index = header_int_index + 1
1959
1960	END SUBROUTINE wrd_mpi_io_int
1961
1962
1963	!--------------------------------------------------------------------------------------------------!
1964	! Description:
1965	! ------------
1966	!> To do: Description missing!
1967	!--------------------------------------------------------------------------------------------------!
1968	SUBROUTINE wrd_mpi_io_real( name, value )
1969
1970	IMPLICIT NONE
1971
1972	CHARACTER(LEN=*), INTENT(IN) :: name !<
1973
1974	REAL(wp), INTENT(IN) :: value !<
1975
1976
1977	IF ( header_real_index == max_nr_real ) THEN
1978	STOP '+++ maximum number of REAL entries in restart file header exceeded'
1979	ENDIF
1980
1981	real_names(header_real_index) = name
1982	real_values(header_real_index) = value
1983	header_real_index = header_real_index + 1
1984
1985	END SUBROUTINE wrd_mpi_io_real
1986
1987
1988
1989	!--------------------------------------------------------------------------------------------------!
1990	! Description:
1991	! ------------
1992	!> Write 2d-REAL array with MPI-IO
1993	!--------------------------------------------------------------------------------------------------!
1994	SUBROUTINE wrd_mpi_io_real_2d( name, data )
1995
1996	IMPLICIT NONE
1997
1998	CHARACTER(LEN=*), INTENT(IN) :: name !<
1999
2000	INTEGER(iwp) :: i !<
2001
2002	#if defined( __parallel )
2003	INTEGER, DIMENSION(rd_status_size) :: status !<
2004	#endif
2005
2006	REAL(wp), INTENT(IN), DIMENSION(nysg:nyng,nxlg:nxrg) :: data !<
2007
2008
2009	IF ( header_array_index == max_nr_arrays ) THEN
2010	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2011	CALL message( 'wrd_mpi_io_real_2d', 'PA0585', 1, 2, 0, 6, 0 )
2012	ENDIF
2013
2014	array_names(header_array_index) = name
2015	array_offset(header_array_index) = array_position
2016	header_array_index = header_array_index + 1
2017
2018	IF ( include_total_domain_boundaries ) THEN
2019	!
2020	!-- Prepare output with outer boundaries
2021	DO i = iog%nxl, iog%nxr
2022	array_2d(i,iog%nys:iog%nyn) = data(iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
2023	ENDDO
2024
2025	ELSE
2026	!
2027	!-- Prepare output without outer boundaries
2028	DO i = nxl,nxr
2029	array_2d(i,iog%nys:iog%nyn) = data(nys:nyn,i)
2030	ENDDO
2031
2032	ENDIF
2033
2034	#if defined( __parallel )
2035	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
2036	IF ( sm_io%iam_io_pe ) THEN
2037	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL, ierr )
2038	CALL MPI_FILE_WRITE_ALL( fh, array_2d, SIZE( array_2d), MPI_REAL, status, ierr )
2039	ENDIF
2040	CALL sm_io%sm_node_barrier()
2041	#else
2042	CALL posix_lseek( fh, array_position )
2043	CALL posix_write( fh, array_2d, SIZE( array_2d ) )
2044	#endif
2045	!
2046	!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
2047	!-- Maybe a compiler problem.
2048	array_position = array_position + ( INT( iog%ny, KIND=rd_offset_kind ) + 1 ) * &
2049	( INT( iog%nx, KIND=rd_offset_kind ) + 1 ) * wp
2050
2051	END SUBROUTINE wrd_mpi_io_real_2d
2052
2053
2054
2055	!--------------------------------------------------------------------------------------------------!
2056	! Description:
2057	! ------------
2058	!> Write 2d-INTEGER array with MPI-IO
2059	!--------------------------------------------------------------------------------------------------!
2060	SUBROUTINE wrd_mpi_io_int_2d( name, data )
2061
2062	IMPLICIT NONE
2063
2064	CHARACTER(LEN=*), INTENT(IN) :: name !<
2065
2066	INTEGER(iwp) :: i !<
2067	INTEGER(iwp) :: j !<
2068
2069	#if defined( __parallel )
2070	INTEGER, DIMENSION(rd_status_size) :: status !<
2071	#endif
2072	INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:,:) :: data !<
2073
2074
2075	IF ( header_array_index == max_nr_arrays ) THEN
2076	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2077	CALL message( 'wrd_mpi_io_int_2d', 'PA0585', 1, 2, 0, 6, 0 )
2078	ENDIF
2079
2080	array_names(header_array_index) = name
2081	array_offset(header_array_index) = array_position
2082	header_array_index = header_array_index + 1
2083
2084	IF ( ( nxr - nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) ) THEN
2085	!
2086	!-- Integer arrays with ghost layers are not implemented yet. These kind of arrays would be
2087	!-- dimensioned in the caller subroutine as
2088	!-- INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) :: data
2089	message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be written to restart ' // &
2090	'file is defined with illegal dimensions in the PALM code'
2091	CALL message( 'wrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
2092
2093	ELSEIF ( ( nxr-nxl+1 ) == SIZE( data, 2 ) ) THEN
2094	!
2095	!-- INTEGER input array without ghost layers.
2096	!-- This kind of array is dimensioned in the caller subroutine as
2097	!-- INTEGER, DIMENSION(nys:nyn,nxl:nxr) :: data
2098	DO j = nys, nyn
2099	DO i = nxl, nxr
2100	array_2di(i,j) = data(j-nys+1,i-nxl+1)
2101	ENDDO
2102	ENDDO
2103	#if defined( __parallel )
2104	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
2105	IF ( sm_io%iam_io_pe ) THEN
2106	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native', &
2107	MPI_INFO_NULL, ierr )
2108	CALL MPI_FILE_WRITE_ALL( fh, array_2di, SIZE( array_2di ), MPI_INTEGER, status, ierr )
2109	ENDIF
2110	CALL sm_io%sm_node_barrier()
2111	#else
2112	CALL posix_lseek( fh, array_position )
2113	CALL posix_write( fh, array_2di, SIZE( array_2di ) )
2114	#endif
2115	!
2116	!-- Type conversion required, otherwise rigth hand side brackets are calculated assuming 4 byte
2117	!-- INT. Maybe a compiler problem.
2118	array_position = array_position + INT( (ny+1), KIND=rd_offset_kind ) * &
2119	INT( (nx+1), KIND=rd_offset_kind ) * 4
2120
2121	ELSE
2122
2123	message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be written to restart ' // &
2124	'file is defined with illegal dimensions in the PALM code'
2125	CALL message( 'wrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
2126
2127	ENDIF
2128
2129	END SUBROUTINE wrd_mpi_io_int_2d
2130
2131
2132
2133	!--------------------------------------------------------------------------------------------------!
2134	! Description:
2135	! ------------
2136	!> Write 3d-INTEGER*4 array with MPI-IO
2137	!--------------------------------------------------------------------------------------------------!
2138	SUBROUTINE wrd_mpi_io_int4_3d( name, data )
2139
2140	IMPLICIT NONE
2141
2142	CHARACTER(LEN=*), INTENT(IN) :: name !<
2143
2144	INTEGER(iwp) :: i !<
2145	#if defined( __parallel )
2146	INTEGER, DIMENSION(rd_status_size) :: status !<
2147	#endif
2148	INTEGER(isp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
2149
2150
2151	IF ( header_array_index == max_nr_arrays ) THEN
2152	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2153	CALL message( 'wrd_mpi_io_int4_3d', 'PA0585', 1, 2, 0, 6, 0 )
2154	ENDIF
2155
2156	array_names(header_array_index) = name
2157	array_offset(header_array_index) = array_position
2158	header_array_index = header_array_index + 1
2159
2160	IF ( include_total_domain_boundaries ) THEN
2161	!
2162	!-- Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
2163	!-- dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
2164	!-- index order of the array in the same way, i.e. the first dimension should be along x and the
2165	!-- second along y. For this reason, the original PALM data need to be swaped.
2166	DO i = iog%nxl, iog%nxr
2167	array_3di4(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
2168	ENDDO
2169
2170	ELSE
2171	!
2172	!-- Prepare output of 3d-REAL-array without ghost layers
2173	DO i = nxl, nxr
2174	array_3di4(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
2175	ENDDO
2176
2177	ENDIF
2178	#if defined( __parallel )
2179	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
2180	IF ( sm_io%iam_io_pe ) THEN
2181	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_3di4, 'native', MPI_INFO_NULL, ierr )
2182	CALL MPI_FILE_WRITE_ALL( fh, array_3di4, SIZE( array_3di4 ), MPI_INTEGER, status, ierr )
2183	ENDIF
2184	CALL sm_io%sm_node_barrier()
2185	#else
2186	CALL posix_lseek( fh, array_position )
2187	CALL posix_write( fh, array_3di4, SIZE( array_3di4 ) )
2188	#endif
2189	!
2190	!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
2191	!-- Maybe a compiler problem.
2192	array_position = array_position + INT( (nz+2), KIND = rd_offset_kind ) * &
2193	INT( (iog%ny+1), KIND = rd_offset_kind ) * &
2194	INT( (iog%nx+1), KIND = rd_offset_kind ) * isp
2195
2196	END SUBROUTINE wrd_mpi_io_int4_3d
2197
2198
2199
2200	!--------------------------------------------------------------------------------------------------!
2201	! Description:
2202	! ------------
2203	!> Write 3d-INTEGER*8 array with MPI-IO
2204	!--------------------------------------------------------------------------------------------------!
2205	SUBROUTINE wrd_mpi_io_int8_3d( name, data )
2206
2207	IMPLICIT NONE
2208
2209	CHARACTER(LEN=*), INTENT(IN) :: name !<
2210
2211	INTEGER(iwp) :: i !<
2212	#if defined( __parallel )
2213	INTEGER, DIMENSION(rd_status_size) :: status !<
2214	#endif
2215	INTEGER(idp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
2216
2217
2218	IF ( header_array_index == max_nr_arrays ) THEN
2219	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2220	CALL message( 'wrd_mpi_io_int8_3d', 'PA0585', 1, 2, 0, 6, 0 )
2221	ENDIF
2222
2223	array_names(header_array_index) = name
2224	array_offset(header_array_index) = array_position
2225	header_array_index = header_array_index + 1
2226
2227	IF ( include_total_domain_boundaries ) THEN
2228	!
2229	!-- Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
2230	!-- dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
2231	!-- index order of the array in the same way, i.e. the first dimension should be along x and the
2232	!-- second along y. For this reason, the original PALM data need to be swaped.
2233	DO i = iog%nxl, iog%nxr
2234	array_3di8(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
2235	ENDDO
2236
2237	ELSE
2238	!
2239	!-- Prepare output of 3d-REAL-array without ghost layers
2240	DO i = nxl, nxr
2241	array_3di8(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
2242	ENDDO
2243
2244	ENDIF
2245	#if defined( __parallel )
2246	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
2247	IF ( sm_io%iam_io_pe ) THEN
2248	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER8, ft_3di8, 'native', MPI_INFO_NULL, ierr )
2249	CALL MPI_FILE_WRITE_ALL( fh, array_3di8, SIZE( array_3di8 ), MPI_INTEGER8, status, ierr )
2250	ENDIF
2251	CALL sm_io%sm_node_barrier()
2252	#else
2253	CALL posix_lseek( fh, array_position )
2254	CALL posix_write( fh, array_3di8, SIZE( array_3di8 ) )
2255	#endif
2256	!
2257	!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
2258	!-- Maybe a compiler problem.
2259	array_position = array_position + INT( (nz+2), KIND = rd_offset_kind ) * &
2260	INT( (iog%ny+1), KIND = rd_offset_kind ) * &
2261	INT( (iog%nx+1), KIND = rd_offset_kind ) * dp
2262
2263	END SUBROUTINE wrd_mpi_io_int8_3d
2264
2265
2266
2267	!--------------------------------------------------------------------------------------------------!
2268	! Description:
2269	! ------------
2270	!> Write 3d-REAL array with MPI-IO
2271	!--------------------------------------------------------------------------------------------------!
2272	SUBROUTINE wrd_mpi_io_real_3d( name, data )
2273
2274	IMPLICIT NONE
2275
2276	CHARACTER(LEN=*), INTENT(IN) :: name !<
2277
2278	INTEGER(iwp) :: i !<
2279	#if defined( __parallel )
2280	INTEGER, DIMENSION(rd_status_size) :: status !<
2281	#endif
2282	REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
2283
2284
2285	IF ( header_array_index == max_nr_arrays ) THEN
2286	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2287	CALL message( 'wrd_mpi_io_real_3d', 'PA0585', 1, 2, 0, 6, 0 )
2288	ENDIF
2289
2290	array_names(header_array_index) = name
2291	array_offset(header_array_index) = array_position
2292	header_array_index = header_array_index + 1
2293
2294	IF ( include_total_domain_boundaries ) THEN
2295	!
2296	!-- Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
2297	!-- dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
2298	!-- index order of the array in the same way, i.e. the first dimension should be along x and the
2299	!-- second along y. For this reason, the original PALM data need to be swaped.
2300	DO i = iog%nxl, iog%nxr
2301	array_3d(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
2302	ENDDO
2303
2304	ELSE
2305	!
2306	!-- Prepare output of 3d-REAL-array without ghost layers
2307	DO i = nxl, nxr
2308	array_3d(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
2309	ENDDO
2310
2311	ENDIF
2312	#if defined( __parallel )
2313	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
2314	IF ( sm_io%iam_io_pe ) THEN
2315	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL, ierr )
2316	CALL MPI_FILE_WRITE_ALL( fh, array_3d, SIZE( array_3d ), MPI_REAL, status, ierr )
2317	ENDIF
2318	CALL sm_io%sm_node_barrier()
2319	#else
2320	CALL posix_lseek( fh, array_position )
2321	CALL posix_write( fh, array_3d, SIZE( array_3d ) )
2322	#endif
2323	!
2324	!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
2325	!-- Maybe a compiler problem.
2326	array_position = array_position + INT( (nz+2), KIND = rd_offset_kind ) * &
2327	INT( (iog%ny+1), KIND = rd_offset_kind ) * &
2328	INT( (iog%nx+1), KIND = rd_offset_kind ) * wp
2329
2330	IF ( debug_output ) THEN
2331	WRITE (9,*) 'array_position real3d ', TRIM( name ), ' ', array_position
2332	FLUSH( 9 )
2333	ENDIF
2334
2335	END SUBROUTINE wrd_mpi_io_real_3d
2336
2337
2338
2339	!--------------------------------------------------------------------------------------------------!
2340	! Description:
2341	! ------------
2342	!> Write 3d-REAL soil array with MPI-IO.
2343	!> nzb_soil, nzt_soil are located in the module land_surface_model_mod. Since Fortran does not allow
2344	!> cross referencing of module variables, it is required to pass these variables as arguments.
2345	!--------------------------------------------------------------------------------------------------!
2346	SUBROUTINE wrd_mpi_io_real_3d_soil( name, data, nzb_soil, nzt_soil )
2347
2348	IMPLICIT NONE
2349
2350	CHARACTER(LEN=*), INTENT(IN) :: name !<
2351
2352	INTEGER(iwp) :: i !<
2353	INTEGER, INTENT(IN) :: nzb_soil !<
2354	INTEGER, INTENT(IN) :: nzt_soil !<
2355
2356	#if defined( __parallel )
2357	INTEGER, DIMENSION(rd_status_size) :: status !<
2358	#endif
2359
2360	REAL(wp), INTENT(IN), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) :: data !<
2361
2362
2363	IF ( header_array_index == max_nr_arrays ) THEN
2364	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2365	CALL message( 'wrd_mpi_io_real_3d_soil', 'PA0585', 1, 2, 0, 6, 0 )
2366	ENDIF
2367
2368	array_names(header_array_index) = name
2369	array_offset(header_array_index) = array_position
2370	header_array_index = header_array_index + 1
2371
2372	CALL rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
2373
2374	IF ( include_total_domain_boundaries) THEN
2375	!
2376	!-- Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
2377	!-- dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
2378	!-- index order of the array in the same way, i.e. the first dimension should be along x and the
2379	!-- second along y. For this reason, the original PALM data need to be swaped.
2380	DO i = iog%nxl, iog%nxr
2381	array_3d_soil(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
2382	ENDDO
2383
2384	ELSE
2385	!
2386	!-- Prepare output of 3d-REAL-array without ghost layers
2387	DO i = nxl, nxr
2388	array_3d_soil(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
2389	ENDDO
2390
2391	ENDIF
2392	#if defined( __parallel )
2393	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
2394	IF ( sm_io%iam_io_pe ) THEN
2395	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native', MPI_INFO_NULL, &
2396	ierr )
2397	CALL MPI_FILE_WRITE_ALL( fh, array_3d_soil, SIZE( array_3d_soil ), MPI_REAL, status, ierr )
2398	ENDIF
2399	CALL sm_io%sm_node_barrier()
2400
2401	IF ( sm_io%is_sm_active() ) THEN
2402	CALL MPI_WIN_FREE( win_3ds, ierr )
2403	ELSE
2404	DEALLOCATE( array_3d_soil )
2405	ENDIF
2406	IF ( sm_io%iam_io_pe ) THEN
2407	CALL MPI_TYPE_FREE( ft_3dsoil, ierr )
2408	ENDIF
2409	#else
2410	CALL posix_lseek( fh, array_position )
2411	CALL posix_write( fh, array_3d_soil, SIZE( array_3d_soil ) )
2412	DEALLOCATE( array_3d_soil )
2413	#endif
2414	!
2415	!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
2416	!-- Maybe a compiler problem.
2417	array_position = array_position + INT( (nzt_soil-nzb_soil+1), KIND = rd_offset_kind ) * &
2418	INT( (iog%ny+1), KIND = rd_offset_kind ) * &
2419	INT( (iog%nx+1), KIND = rd_offset_kind ) * wp
2420
2421	END SUBROUTINE wrd_mpi_io_real_3d_soil
2422
2423
2424
2425	!--------------------------------------------------------------------------------------------------!
2426	! Description:
2427	! ------------
2428	!> Write CHARATCTER with MPI-IO
2429	!--------------------------------------------------------------------------------------------------!
2430	SUBROUTINE wrd_mpi_io_char( name, text )
2431
2432	IMPLICIT NONE
2433
2434	CHARACTER(LEN=128) :: lo_line !<
2435	CHARACTER(LEN=*), INTENT(IN) :: name !<
2436	CHARACTER(LEN=*), INTENT(IN) :: text !<
2437
2438
2439	IF ( header_char_index == max_nr_char ) THEN
2440	STOP '+++ maximum number of CHARACTER entries in restart file header exceeded'
2441	ENDIF
2442
2443	lo_line = name
2444	lo_line(33:) = text
2445	text_lines(header_char_index) = lo_line
2446	header_char_index = header_char_index + 1
2447
2448	END SUBROUTINE wrd_mpi_io_char
2449
2450
2451
2452	!--------------------------------------------------------------------------------------------------!
2453	! Description:
2454	! ------------
2455	!> Write LOGICAL with MPI-IO
2456	!--------------------------------------------------------------------------------------------------!
2457	SUBROUTINE wrd_mpi_io_logical( name, value )
2458
2459	IMPLICIT NONE
2460
2461	CHARACTER(LEN=*), INTENT(IN) :: name !<
2462
2463	INTEGER(iwp) :: logical_as_integer !<
2464
2465	LOGICAL, INTENT(IN) :: value !<
2466
2467
2468	IF ( value ) THEN
2469	logical_as_integer = 1
2470	ELSE
2471	logical_as_integer = 0
2472	ENDIF
2473
2474	CALL wrd_mpi_io_int( name, logical_as_integer )
2475
2476	END SUBROUTINE wrd_mpi_io_logical
2477
2478
2479
2480	!--------------------------------------------------------------------------------------------------!
2481	! Description:
2482	! ------------
2483	!> Write 3D-LOGICAL with MPI-IO
2484	!--------------------------------------------------------------------------------------------------!
2485	SUBROUTINE wrd_mpi_io_logical_3d( name, data )
2486
2487	IMPLICIT NONE
2488
2489	CHARACTER(LEN=*), INTENT(IN) :: name !<
2490
2491	INTEGER(iwp) :: i !< loop index
2492	INTEGER(iwp) :: j !< loop index
2493	INTEGER(iwp) :: k !< loop index
2494
2495	INTEGER(iwp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: logical_as_integer !<
2496
2497	LOGICAL, INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: data !<
2498
2499
2500	DO i = nxlg, nxrg
2501	DO j = nysg, nyng
2502	DO k = nzb, nzt+1
2503	IF ( data(k,j,i) ) THEN
2504	logical_as_integer(k,j,i) = 1
2505	ELSE
2506	logical_as_integer(k,j,i) = 0
2507	ENDIF
2508	ENDDO
2509	ENDDO
2510	ENDDO
2511
2512	CALL wrd_mpi_io_int4_3d( name, logical_as_integer )
2513
2514	END SUBROUTINE wrd_mpi_io_logical_3d
2515
2516
2517
2518	!--------------------------------------------------------------------------------------------------!
2519	! Description:
2520	! ------------
2521	!> Read 1d-REAL global array with MPI-IO.
2522	!> Array contains identical data on all PEs.
2523	!--------------------------------------------------------------------------------------------------!
2524	SUBROUTINE rrd_mpi_io_global_array_real_1d( name, data )
2525
2526	IMPLICIT NONE
2527
2528	CHARACTER(LEN=*), INTENT(IN) :: name !<
2529
2530	INTEGER(iwp) :: i !<
2531	INTEGER(KIND=rd_offset_kind) :: offset !<
2532
2533	#if defined( __parallel )
2534	INTEGER, DIMENSION(rd_status_size) :: status !<
2535	#endif
2536
2537	LOGICAL :: found !<
2538
2539	REAL(KIND=wp), INTENT(INOUT), DIMENSION(:) :: data !<
2540
2541
2542	offset = 0
2543	found = .FALSE.
2544
2545	DO i = 1, tgh%nr_arrays
2546	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
2547	array_position = array_offset(i)
2548	found = .TRUE.
2549	EXIT
2550	ENDIF
2551	ENDDO
2552
2553
2554	IF ( found ) THEN
2555
2556	!
2557	!-- Set default view
2558	#if defined( __parallel )
2559	IF ( cyclic_fill_mode ) THEN !kk This may be the general solution for all cases
2560	IF ( pe_active_for_read ) THEN
2561	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2562	CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
2563	CALL MPI_FILE_READ_ALL( fh, data, SIZE( data ), MPI_REAL, status, ierr )
2564	ENDIF
2565	CALL MPI_BCAST( data, SIZE( data ), MPI_REAL, 0, comm2d, ierr )
2566	ELSE
2567	IF( sm_io%iam_io_pe ) THEN
2568	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2569	ENDIF
2570	IF ( myid == 0 ) THEN
2571	CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
2572	CALL MPI_FILE_READ( fh, data, SIZE( data ), MPI_REAL, status, ierr )
2573	ENDIF
2574	CALL MPI_BCAST( data, SIZE( data ), MPI_REAL, 0, comm2d, ierr )
2575	ENDIF
2576	#else
2577	CALL posix_lseek( fh, array_position )
2578	CALL posix_read( fh, data, SIZE( data ) )
2579	#endif
2580
2581	ELSE
2582
2583	message_string = '1d/2d/3d/4d-REAL global array "' // TRIM( name ) // '" not found in ' // &
2584	'restart file'
2585	CALL message( 'rrd_mpi_io_global_array_real_1d', 'PA0722', 3, 2, 0, 6, 0 )
2586
2587	ENDIF
2588
2589	END SUBROUTINE rrd_mpi_io_global_array_real_1d
2590
2591
2592
2593	!--------------------------------------------------------------------------------------------------!
2594	! Description:
2595	! ------------
2596	!> Read 2d-REAL global array with MPI-IO.
2597	!> Array contains identical data on all PEs.
2598	!--------------------------------------------------------------------------------------------------!
2599	SUBROUTINE rrd_mpi_io_global_array_real_2d( name, data )
2600
2601	IMPLICIT NONE
2602
2603	CHARACTER(LEN=*), INTENT(IN) :: name !<
2604
2605	INTEGER, DIMENSION(1) :: bufshape !<
2606
2607	REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET :: data !<
2608	REAL(KIND=wp), POINTER, DIMENSION(:) :: buf !<
2609
2610	TYPE(C_PTR) :: c_data !<
2611
2612
2613	c_data = C_LOC( data )
2614	bufshape(1) = SIZE( data )
2615	CALL C_F_POINTER( c_data, buf, bufshape )
2616
2617	CALL rrd_mpi_io_global_array_real_1d( name, buf )
2618
2619	END SUBROUTINE rrd_mpi_io_global_array_real_2d
2620
2621
2622
2623	!--------------------------------------------------------------------------------------------------!
2624	! Description:
2625	! ------------
2626	!> Read 3d-REAL global array with MPI-IO.
2627	!> Array contains identical data on all PEs.
2628	!--------------------------------------------------------------------------------------------------!
2629	SUBROUTINE rrd_mpi_io_global_array_real_3d( name, data )
2630
2631	IMPLICIT NONE
2632
2633	CHARACTER(LEN=*), INTENT(IN) :: name !<
2634
2635	INTEGER, DIMENSION(1) :: bufshape !<
2636
2637	REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET :: data !<
2638	REAL(KIND=wp), POINTER, DIMENSION(:) :: buf !<
2639
2640	TYPE(C_PTR) :: c_data !<
2641
2642
2643	c_data = C_LOC( data )
2644	bufshape(1) = SIZE( data )
2645	CALL C_F_POINTER( c_data, buf, bufshape )
2646
2647	CALL rrd_mpi_io_global_array_real_1d( name, buf )
2648
2649	END SUBROUTINE rrd_mpi_io_global_array_real_3d
2650
2651
2652
2653	!--------------------------------------------------------------------------------------------------!
2654	! Description:
2655	! ------------
2656	!> Read 4d-REAL global array with MPI-IO.
2657	!> Array contains identical data on all PEs.
2658	!--------------------------------------------------------------------------------------------------!
2659	SUBROUTINE rrd_mpi_io_global_array_real_4d( name, data )
2660
2661	IMPLICIT NONE
2662
2663	CHARACTER(LEN=*), INTENT(IN) :: name !<
2664
2665	INTEGER, DIMENSION(1) :: bufshape !<
2666
2667	REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET :: data !<
2668	REAL(KIND=wp), POINTER, DIMENSION(:) :: buf !<
2669
2670	TYPE(C_PTR) :: c_data !<
2671
2672
2673	c_data = C_LOC( data )
2674	bufshape(1) = SIZE( data)
2675	CALL C_F_POINTER( c_data, buf, bufshape )
2676
2677	CALL rrd_mpi_io_global_array_real_1d( name, buf )
2678
2679	END SUBROUTINE rrd_mpi_io_global_array_real_4d
2680
2681
2682
2683	!--------------------------------------------------------------------------------------------------!
2684	! Description:
2685	! ------------
2686	!> Read 1d-INTEGER global array with MPI-IO.
2687	!> Array contains identical data on all PEs.
2688	!--------------------------------------------------------------------------------------------------!
2689	SUBROUTINE rrd_mpi_io_global_array_int_1d( name, data )
2690
2691	IMPLICIT NONE
2692
2693	CHARACTER(LEN=*), INTENT(IN) :: name !<
2694
2695	INTEGER(iwp) :: i !<
2696	INTEGER(KIND=rd_offset_kind) :: offset !<
2697
2698	#if defined( __parallel )
2699	INTEGER, DIMENSION(rd_status_size) :: status !<
2700	#endif
2701	INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:) :: data !<
2702
2703	LOGICAL :: found !<
2704
2705
2706	offset = 0
2707	found = .FALSE.
2708
2709	DO i = 1, tgh%nr_arrays
2710	IF ( TRIM( array_names(i) ) == TRIM( name ) ) THEN
2711	array_position = array_offset(i)
2712	found = .TRUE.
2713	EXIT
2714	ENDIF
2715	ENDDO
2716
2717	IF ( found ) THEN
2718	!
2719	!-- Set default view
2720	#if defined( __parallel )
2721	IF ( cyclic_fill_mode ) THEN !kk This may be the general solution for all cases
2722	IF ( pe_active_for_read ) THEN
2723	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2724	CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
2725	CALL MPI_FILE_READ_ALL( fh, data, SIZE( data), MPI_INTEGER, status, ierr )
2726	ENDIF
2727	CALL MPI_BCAST( data, SIZE( data ), MPI_INTEGER, 0, comm2d, ierr )
2728	ELSE
2729	IF( sm_io%iam_io_pe ) THEN
2730	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2731	ENDIF
2732	IF ( myid == 0 ) THEN
2733	CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
2734	CALL MPI_FILE_READ( fh, data, SIZE( data), MPI_INTEGER, status, ierr )
2735	ENDIF
2736	CALL MPI_BCAST( data, SIZE( data ), MPI_INTEGER, 0, comm2d, ierr )
2737	ENDIF
2738	#else
2739	CALL posix_lseek( fh, array_position )
2740	CALL posix_read( fh, data, SIZE( data ) )
2741	#endif
2742	ELSE
2743
2744	message_string = '1d-INTEGER global array "' // TRIM( name ) // '" not found in ' // &
2745	'restart file'
2746	CALL message( 'rrd_mpi_io_global_array_int_1d', 'PA0722', 3, 2, 0, 6, 0 )
2747
2748	ENDIF
2749
2750	END SUBROUTINE rrd_mpi_io_global_array_int_1d
2751
2752
2753
2754	!--------------------------------------------------------------------------------------------------!
2755	! Description:
2756	! ------------
2757	!> Write 1d-REAL global array with MPI-IO.
2758	!> Array contains identical data on all PEs.
2759	!--------------------------------------------------------------------------------------------------!
2760	SUBROUTINE wrd_mpi_io_global_array_real_1d( name, data )
2761
2762	IMPLICIT NONE
2763
2764	CHARACTER(LEN=*), INTENT(IN) :: name !<
2765
2766	INTEGER(KIND=rd_offset_kind) :: offset !<
2767
2768	#if defined( __parallel )
2769	INTEGER, DIMENSION(rd_status_size) :: status !<
2770	#endif
2771
2772	REAL(KIND=wp), INTENT(IN), DIMENSION(:) :: data !<
2773
2774
2775	offset = 0
2776
2777	IF ( header_array_index == max_nr_arrays ) THEN
2778	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2779	CALL message( 'wrd_mpi_io_global_array_real_1d', 'PA0585', 1, 2, 0, 6, 0 )
2780	ENDIF
2781
2782	array_names(header_array_index) = name
2783	array_offset(header_array_index) = array_position
2784	header_array_index = header_array_index + 1
2785
2786	!
2787	!-- Set default view
2788	#if defined( __parallel )
2789	IF ( sm_io%iam_io_pe ) THEN
2790	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2791	ENDIF
2792	!
2793	!-- Only PE 0 writes replicated data
2794	IF ( myid == 0 ) THEN
2795	CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
2796	CALL MPI_FILE_WRITE( fh, data, SIZE( data ), MPI_REAL, status, ierr )
2797	ENDIF
2798	#else
2799	CALL posix_lseek( fh, array_position )
2800	CALL posix_write( fh, data, SIZE( data ) )
2801	#endif
2802	array_position = array_position + SIZE( data ) * wp
2803
2804	END SUBROUTINE wrd_mpi_io_global_array_real_1d
2805
2806
2807
2808	!--------------------------------------------------------------------------------------------------!
2809	! Description:
2810	! ------------
2811	!> Write 2d-REAL global array with MPI-IO.
2812	!> Array contains identical data on all PEs.
2813	!--------------------------------------------------------------------------------------------------!
2814	SUBROUTINE wrd_mpi_io_global_array_real_2d( name, data )
2815
2816	IMPLICIT NONE
2817
2818	CHARACTER(LEN=*), INTENT(IN) :: name !<
2819
2820	INTEGER, DIMENSION(1) :: bufshape !<
2821
2822	REAL(KIND=wp), POINTER, DIMENSION(:) :: buf !<
2823	REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET :: data !<
2824
2825	TYPE(C_PTR) :: c_data !<
2826
2827
2828	c_data = C_LOC( data )
2829	bufshape(1) = SIZE( data)
2830	CALL C_F_POINTER( c_data, buf, bufshape )
2831
2832	CALL wrd_mpi_io_global_array_real_1d( name, buf )
2833
2834	END SUBROUTINE wrd_mpi_io_global_array_real_2d
2835
2836
2837
2838	!--------------------------------------------------------------------------------------------------!
2839	! Description:
2840	! ------------
2841	!> Write 3d-REAL global array with MPI-IO.
2842	!> Array contains identical data on all PEs.
2843	!--------------------------------------------------------------------------------------------------!
2844	SUBROUTINE wrd_mpi_io_global_array_real_3d( name, data )
2845
2846	IMPLICIT NONE
2847
2848	CHARACTER(LEN=*), INTENT(IN) :: name !<
2849
2850	INTEGER, DIMENSION(1) :: bufshape !<
2851
2852	REAL(KIND=wp), POINTER, DIMENSION(:) :: buf !<
2853	REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET :: data !<
2854
2855	TYPE(C_PTR) :: c_data !<
2856
2857
2858	c_data = C_LOC( data )
2859	bufshape(1) = SIZE( data )
2860	CALL C_F_POINTER( c_data, buf, bufshape )
2861
2862	CALL wrd_mpi_io_global_array_real_1d( name, buf )
2863
2864	END SUBROUTINE wrd_mpi_io_global_array_real_3d
2865
2866
2867
2868	!--------------------------------------------------------------------------------------------------!
2869	! Description:
2870	! ------------
2871	!> Write 4d-REAL global array with MPI-IO.
2872	!> Array contains identical data on all PEs.
2873	!--------------------------------------------------------------------------------------------------!
2874	SUBROUTINE wrd_mpi_io_global_array_real_4d( name, data )
2875
2876	IMPLICIT NONE
2877
2878	CHARACTER(LEN=*), INTENT(IN) :: name !<
2879
2880	INTEGER, DIMENSION(1) :: bufshape !<
2881
2882	REAL(KIND=wp), POINTER, DIMENSION(:) :: buf !<
2883	REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET :: data !<
2884
2885	TYPE(C_PTR) :: c_data !<
2886
2887
2888	c_data = C_LOC( data )
2889	bufshape(1) = SIZE( data)
2890	CALL C_F_POINTER( c_data, buf, bufshape )
2891
2892	CALL wrd_mpi_io_global_array_real_1d( name, buf )
2893
2894	END SUBROUTINE wrd_mpi_io_global_array_real_4d
2895
2896
2897
2898	!--------------------------------------------------------------------------------------------------!
2899	! Description:
2900	! ------------
2901	!> Write 1d-INTEGER global array with MPI-IO.
2902	!> Array contains identical data on all PEs.
2903	!--------------------------------------------------------------------------------------------------!
2904	SUBROUTINE wrd_mpi_io_global_array_int_1d( name, data )
2905
2906	IMPLICIT NONE
2907
2908	CHARACTER(LEN=*), INTENT(IN) :: name !<
2909
2910	INTEGER(KIND=rd_offset_kind) :: offset !<
2911
2912	INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:) :: data !<
2913	#if defined( __parallel )
2914	INTEGER, DIMENSION(rd_status_size) :: status !<
2915	#endif
2916
2917	IF ( header_array_index == max_nr_arrays ) THEN
2918	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
2919	CALL message( 'wrd_mpi_io_global_array_int_1d', 'PA0585', 1, 2, 0, 6, 0 )
2920	ENDIF
2921
2922	offset = 0
2923	array_names(header_array_index) = name
2924	array_offset(header_array_index) = array_position
2925	header_array_index = header_array_index + 1
2926
2927	!
2928	!-- Set default view
2929	#if defined( __parallel )
2930	IF ( sm_io%iam_io_pe ) THEN
2931	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2932	ENDIF
2933	!
2934	!-- Only PE 0 writes replicated data
2935	IF ( myid == 0 ) THEN
2936	CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
2937	CALL MPI_FILE_WRITE( fh, data, SIZE( data), MPI_INTEGER, status, ierr )
2938	ENDIF
2939	#else
2940	CALL posix_lseek( fh, array_position )
2941	CALL posix_write( fh, data, SIZE( data ) )
2942	#endif
2943	array_position = array_position + SIZE( data ) * 4
2944
2945	END SUBROUTINE wrd_mpi_io_global_array_int_1d
2946
2947
2948
2949	!--------------------------------------------------------------------------------------------------!
2950	! Description:
2951	! ------------
2952	!> Read particle data array with MPI-IO.
2953	!--------------------------------------------------------------------------------------------------!
2954	SUBROUTINE rrd_mpi_io_particles( name, prt_global_index )
2955
2956	IMPLICIT NONE
2957
2958	CHARACTER(LEN=*), INTENT(IN) :: name !<
2959	INTEGER(idp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: prt_global_index !<
2960
2961	INTEGER(iwp) :: array_size !<
2962	INTEGER(iwp) :: byte_column !<
2963	INTEGER(iwp) :: i !<
2964	INTEGER(iwp) :: ind !<
2965	INTEGER(iwp) :: j !<
2966	INTEGER(iwp) :: k !<
2967	INTEGER(iwp) :: n !<
2968	INTEGER(iwp) :: particle_size !<
2969
2970	INTEGER(KIND=rd_offset_kind) :: disp !<
2971	INTEGER(KIND=rd_offset_kind) :: offset !<
2972	INTEGER(KIND=rd_offset_kind) :: prt_nr_bytes !<
2973
2974	LOGICAL :: found !<
2975
2976	REAL(dp) :: rr !< there is no data type INTEGER*8 in MPI
2977	REAL(dp) :: rs !< use REAL*8 to compute max offset
2978
2979	TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET :: prt_data !<
2980
2981	#if defined( __parallel )
2982	INTEGER, DIMENSION(rd_status_size) :: status !<
2983	#else
2984	TYPE(C_PTR) :: buf
2985	#endif
2986
2987	found = .FALSE.
2988
2989	DO i = 1, tgh%nr_arrays
2990	IF ( TRIM(array_names(i)) == TRIM( name ) ) THEN
2991	array_position = array_offset(i)
2992	found = .TRUE.
2993	EXIT
2994	ENDIF
2995	ENDDO
2996
2997	IF ( found ) THEN
2998
2999	offset = 0
3000
3001	particle_size = STORAGE_SIZE(zero_particle) / 8 ! 8 here means number of bits/byte, NOT wp
3002
3003	array_size = 0
3004	DO i = nxl, nxr
3005	DO j = nys, nyn
3006	array_size = MAX( array_size, SUM(prt_count(:,j,i)) )
3007	ENDDO
3008	ENDDO
3009
3010	ALLOCATE( prt_data(MAX(array_size,1)) )
3011
3012	!
3013	!-- Write columns of particle
3014	#if defined( __parallel )
3015	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
3016	#endif
3017	prt_nr_bytes = 0
3018	DO i = nxl, nxr
3019	DO j = nys, nyn
3020	disp = array_position + prt_global_index(nzb,j,i) * particle_size
3021	byte_column = SUM( prt_count(:,j,i) ) * particle_size
3022	prt_nr_bytes = MAX( disp+byte_column, prt_nr_bytes )
3023
3024	#if defined( __parallel )
3025	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
3026	IF ( byte_column > 0 ) THEN
3027	CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
3028	CALL MPI_FILE_READ( fh, prt_data, byte_column, MPI_BYTE, status, ierr )
3029	ENDIF
3030	CALL sm_io%sm_node_barrier()
3031	#else
3032	buf = C_LOC(prt_data) ! use C_PTR to avaid another overlay in posix interface
3033	CALL posix_lseek( fh, disp )
3034	CALL posix_read( fh, buf, byte_column )
3035	#endif
3036	ind = 1
3037	DO k = nzb, nzt+1
3038	DO n = 1, prt_count(k,j,i)
3039	grid_particles(k,j,i)%particles(n) = prt_data(ind)
3040	ind = ind+1
3041	ENDDO
3042	ENDDO
3043
3044	ENDDO
3045	ENDDO
3046
3047	#if defined( __parallel )
3048	rs = prt_nr_bytes
3049	CALL MPI_ALLREDUCE( rs, rr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm2d, ierr )
3050	prt_nr_bytes = rr
3051	#else
3052	rr = rs
3053	#endif
3054	array_position = prt_nr_bytes
3055
3056	DEALLOCATE( prt_data )
3057
3058	ENDIF
3059
3060	END SUBROUTINE rrd_mpi_io_particles
3061
3062
3063
3064	!--------------------------------------------------------------------------------------------------!
3065	! Description:
3066	! ------------
3067	!> Read 1d-REAL surface data array with MPI-IO.
3068	!> This is a recursive subroutine. In case of cyclic fill mode it may call itself for reading parts
3069	!> of the prerun grid.
3070	!--------------------------------------------------------------------------------------------------!
3071	RECURSIVE SUBROUTINE rrd_mpi_io_surface( name, data, first_index )
3072
3073	IMPLICIT NONE
3074
3075	CHARACTER(LEN=*), INTENT(IN) :: name !<
3076
3077	INTEGER(iwp), OPTIONAL :: first_index !<
3078	INTEGER(iwp) :: i !<
3079	INTEGER(iwp) :: j !<
3080	INTEGER(iwp) :: lo_first_index !<
3081
3082	#if defined( __parallel )
3083	INTEGER(iwp) :: buf_start !<
3084	INTEGER(KIND=rd_offset_kind) :: disp !< displacement of actual indices
3085	INTEGER(iwp) :: ie !<
3086	INTEGER(iwp) :: ind_gb !<
3087	INTEGER(iwp) :: ind_out !<
3088	INTEGER(iwp) :: is !<
3089	INTEGER(iwp) :: n !<
3090	INTEGER(iwp) :: n_trans !<
3091
3092	INTEGER(iwp),DIMENSION(0:numprocs-1) :: lo_index !<
3093	INTEGER, DIMENSION(rd_status_size) :: status !<
3094	#endif
3095	LOGICAL :: found !<
3096
3097	REAL(wp), INTENT(OUT), DIMENSION(:), TARGET :: data !<
3098	#if defined( __parallel )
3099	REAL(wp),DIMENSION(:),ALLOCATABLE :: put_buffer !<
3100	#endif
3101
3102
3103	found = .FALSE.
3104	lo_first_index = 1
3105
3106	IF ( PRESENT( first_index ) ) THEN
3107	lo_first_index = first_index
3108	ENDIF
3109
3110	DO i = 1, tgh%nr_arrays
3111	IF ( TRIM( array_names(i) ) == TRIM( name ) ) THEN
3112	!
3113	!-- ATTENTION: The total_number_of_surface_elements and wp MUST be INTERGER(8).
3114	!-- The compiler (at least Intel) first computes total_number_of_surface_elements*wp
3115	!-- and then does the conversion to INTEGER(8).
3116	!-- This may lead to wrong results when total_number_of_surface_elementswp is > 210**6
3117	array_position = array_offset(i) + ( lo_first_index - 1 ) * &
3118	INT( total_number_of_surface_elements, idp ) * INT( wp, idp )
3119	found = .TRUE.
3120	EXIT
3121	ENDIF
3122	ENDDO
3123
3124	!
3125	!-- In case of 2d-data, name is written only once
3126	IF ( lo_first_index == 1 ) THEN
3127
3128	array_names(header_array_index) = name
3129	array_offset(header_array_index) = array_position
3130	header_array_index = header_array_index + 1
3131
3132	ENDIF
3133
3134	IF ( found ) THEN
3135	IF ( cyclic_fill_mode ) THEN
3136
3137	CALL rrd_mpi_io_surface_cyclic_fill
3138	RETURN
3139
3140	ELSE
3141	#if defined( __parallel )
3142	!
3143	!-- Read data from restart file
3144	CALL sm_io%sm_node_barrier() ! has no effect if I/O on limited number of cores is inactive
3145	IF ( sm_io%iam_io_pe ) THEN
3146	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_surf, 'native', &
3147	MPI_INFO_NULL, ierr )
3148	CALL MPI_FILE_READ_ALL ( fh, array_out, SIZE(array_out), MPI_REAL, status, ierr )
3149	ENDIF
3150	CALL sm_io%sm_node_barrier()
3151
3152	!
3153	!-- Copy data into transfer buffer. Data is organized in a way that only one MPI_PUT to the
3154	!-- respective PE ist required.
3155	ALLOCATE( put_buffer(SUM( transfer_index(4,:) )) )
3156
3157	ind_gb = 1
3158	DO i = 1, SIZE( local_indices, 2 )
3159	ind_out = local_indices(1,i)
3160	DO j = 1, local_indices(2,i)
3161	put_buffer(ind_gb) = array_out(ind_out)
3162	ind_out = ind_out + 1
3163	ind_gb = ind_gb + 1
3164	ENDDO
3165	ENDDO
3166	!
3167	!-- Transfer data from I/O PEs to the respective PEs to which they belong.
3168	CALL MPI_WIN_FENCE( 0, win_surf, ierr )
3169
3170	buf_start = 1
3171	DO n = 0, numprocs-1
3172	n_trans = transfer_index(4,n)
3173	IF ( n_trans > 0 ) THEN
3174	disp = transfer_index(3,n) - 1
3175	CALL MPI_PUT( put_buffer(buf_start), n_trans, MPI_REAL, n, disp, n_trans, MPI_REAL,&
3176	win_surf, ierr)
3177	buf_start = buf_start + n_trans
3178	ENDIF
3179	ENDDO
3180
3181	CALL MPI_WIN_FENCE( 0, win_surf, ierr )
3182	DEALLOCATE( put_buffer )
3183	!
3184	!-- Copy from RMA window into output array (data) to allow transfering data to target PEs.
3185	!-- Check, if the number of surface elements per horizontal gridbox match the index setup.
3186	lo_index = nr_surfaces_in_tb
3187	DO i = nxl, nxr
3188	DO j = nys, nyn
3189	is = lo_index(s_index_in_window(j,i)) + 1
3190	ie = is + m_end_index(j,i) - m_start_index(j,i)
3191	data(m_start_index(j,i):m_end_index(j,i)) = array_1d(is:ie)
3192	lo_index(s_index_in_window(j,i)) = lo_index(s_index_in_window(j,i)) + &
3193	e_end_index(j,i) - e_start_index(j,i) + 1
3194	!
3195	!-- TODO: Test can be removed later.
3196	IF ( e_end_index(j,i)-e_start_index(j,i)+1 /= NINT( array_1d(is-1) ) ) THEN
3197	WRITE( 9, '(A,6I8)' ) 'Nr surface elements does not match ', j, i, &
3198	e_start_index(j,i), e_end_index(j,i), &
3199	e_end_index(j,i)-e_start_index(j,i)+1 , &
3200	NINT( array_1d(is-1) )
3201	FLUSH( 9 )
3202	CALL MPI_ABORT( comm2d, 1, ierr )
3203	ENDIF
3204	ENDDO
3205	ENDDO
3206
3207
3208	#else
3209	CALL posix_lseek( fh, array_position )
3210	CALL posix_read( fh, array_out, SIZE(array_out) )
3211
3212	DO i = nxl, nxr
3213	DO j = nys, nyn
3214	data(m_start_index(j,i):m_end_index(j,i)) = &
3215	array_out(e_start_index(j,i)+1:e_end_index(j,i))
3216	!
3217	!-- TODO: Test can be removed later.
3218	IF ( e_end_index(j,i)-e_start_index(j,i)+1 /= NINT(array_out(e_start_index(j,i))) )&
3219	THEN
3220	WRITE( 9, '(A,6I8)' ) 'Nr surface elements does not match ', j, i, &
3221	e_start_index(j,i), e_end_index(j,i), &
3222	e_end_index(j,i)-e_start_index(j,i)+1, &
3223	NINT( array_out(e_start_index(j,i)) )
3224	FLUSH( 9 )
3225	CALL ABORT()
3226	ENDIF
3227	ENDDO
3228	ENDDO
3229	#endif
3230	ENDIF
3231
3232	ENDIF
3233
3234	CONTAINS
3235
3236	SUBROUTINE rrd_mpi_io_surface_cyclic_fill
3237
3238	IMPLICIT NONE
3239
3240	INTEGER(iwp) :: i !<
3241	INTEGER(iwp) :: ie !<
3242	INTEGER(iwp) :: is !<
3243	INTEGER(iwp) :: i_remote !<
3244	INTEGER(iwp) :: j !<
3245	INTEGER(iwp) :: je !<
3246	INTEGER(iwp) :: js !<
3247	INTEGER(iwp) :: j_remote !<
3248	INTEGER(iwp) :: nval !<
3249	INTEGER(iwp) :: rem_pe !<
3250
3251	#if defined( __parallel )
3252	INTEGER(KIND=MPI_ADDRESS_KIND) :: rem_offs !<
3253	#else
3254	INTEGER(idp) :: rem_offs !<
3255	#endif
3256
3257	REAL(wp), DIMENSION(:), ALLOCATABLE :: c_data !<
3258
3259
3260	!
3261	!-- ATTENTION: This version allows only 1 surface element per grid cell.
3262	!
3263	!-- Activate grid of the smaller prerun, i.e. grid variables like nxl, nxr, nys, nyn and other
3264	!-- values are set according to the prerun settings.
3265	CALL prerun_grid%activate_grid_from_this_class()
3266
3267	IF ( pe_active_for_read ) THEN
3268
3269	IF ( MAXVAL( m_end_index ) <= 0 ) THEN
3270	CALL mainrun_grid%activate_grid_from_this_class()
3271	IF ( debug_output ) THEN
3272	CALL debug_message( 'PE inactive for reading restart or prerun data', 'start' )
3273	ENDIF
3274	RETURN
3275	ENDIF
3276
3277	ALLOCATE( c_data(MAXVAL( m_end_index )) )
3278
3279	!
3280	!-- Recursive CALL of rrd_mpi_io_surface.
3281	!-- rrd_mpi_io_surface is called with cyclic_fill_mode = .FALSE. on the smaller prerun grid.
3282	cyclic_fill_mode = .FALSE.
3283	CALL rrd_mpi_io_surface( name, c_data )
3284	cyclic_fill_mode = .TRUE.
3285
3286	DO i = nxl, nxr
3287	DO j = nys, nyn
3288	rmabuf_2d(j,i) = c_data(c_start_index(j,i))
3289	ENDDO
3290	ENDDO
3291
3292	ENDIF
3293	!
3294	!-- Activate grid of the mainrun, i.e. grid variables like nxl, nxr, nys, nyn and other values
3295	!-- are set according to the mainrun settings.
3296	CALL mainrun_grid%activate_grid_from_this_class()
3297
3298	#if defined( __parallel )
3299	!
3300	!-- Close RMA window to allow remote access
3301	CALL MPI_WIN_FENCE( 0, rmawin_2d, ierr )
3302	#endif
3303
3304	!
3305	!-- After reading surface data on the small grid, map these data in a cyclic way to all respective
3306	!-- grid points of the main run.
3307	is = nxl
3308	ie = nxr
3309	js = nys
3310	je = nyn
3311
3312	DO i = is, ie
3313	DO j = js, je
3314	i_remote = MOD( i, nx_on_file+1 )
3315	j_remote = MOD( j, ny_on_file+1 )
3316	rem_pe = remote_pe(i_remote,j_remote)
3317	rem_offs = rma_offset(i_remote,j_remote)
3318	nval = 1
3319
3320	#if defined( __parallel )
3321	IF ( rem_pe /= myid ) THEN
3322	CALL MPI_GET( data(o_start_index(j,i)), nval, MPI_REAL, rem_pe, rem_offs, nval, &
3323	MPI_REAL, rmawin_2d, ierr)
3324	ELSE
3325	data(o_start_index(j,i)) = rmabuf_2d(j_remote,i_remote)
3326	ENDIF
3327	#else
3328	data(o_start_index(j,i)) = array_2d(i_remote,j_remote)
3329	#endif
3330	ENDDO
3331	ENDDO
3332
3333	#if defined( __parallel )
3334	!
3335	!-- Reopen RMA window to allow filling
3336	CALL MPI_WIN_FENCE( 0, rmawin_2d, ierr )
3337	#endif
3338
3339	IF ( ALLOCATED( c_data ) ) DEALLOCATE( c_data )
3340
3341	END SUBROUTINE rrd_mpi_io_surface_cyclic_fill
3342
3343	END SUBROUTINE rrd_mpi_io_surface
3344
3345
3346
3347	!--------------------------------------------------------------------------------------------------!
3348	! Description:
3349	! ------------
3350	!> Read 2d-REAL surface data array with MPI-IO.
3351	!> These consist of multiple 1d-REAL surface data arrays.
3352	!--------------------------------------------------------------------------------------------------!
3353	SUBROUTINE rrd_mpi_io_surface_2d( name, data )
3354
3355	IMPLICIT NONE
3356
3357	CHARACTER(LEN=*), INTENT(IN) :: name !<
3358
3359	INTEGER(iwp) :: i !<
3360
3361	REAL(wp), INTENT(OUT), DIMENSION(:,:) :: data !<
3362	REAL(wp), DIMENSION(SIZE( data,2)) :: tmp !<
3363
3364
3365	DO i = 1, SIZE( data, 1 )
3366	CALL rrd_mpi_io_surface( name, tmp, first_index = i )
3367	data(i,:) = tmp
3368	ENDDO
3369
3370	!
3371	!-- Comment from Klaus Ketelsen (September 2018)
3372	!-- The intention of the following loop was to let the compiler do the copying on return.
3373	!-- In my understanding it is standard conform to pass the second dimension to a 1d-array inside a
3374	!-- subroutine and the compiler is responsible to generate code for copying. Acually this works fine
3375	!-- for INTENT(IN) variables (wrd_mpi_io_surface_2d). For INTENT(OUT) like in this case the code
3376	!-- works on pgi compiler. But both, the Intel 16 and the Cray compiler show wrong answers using
3377	!-- this loop. That is the reason why the above auxiliary array tmp was introduced.
3378	! DO i = 1, SIZE( data,1)
3379	! CALL rrd_mpi_io_surface( name, data(i,:), first_index = i )
3380	! ENDDO
3381
3382	END SUBROUTINE rrd_mpi_io_surface_2d
3383
3384
3385
3386	!--------------------------------------------------------------------------------------------------!
3387	! Description:
3388	! ------------
3389	!> Write particle data with MPI-IO.
3390	!--------------------------------------------------------------------------------------------------!
3391	SUBROUTINE wrd_mpi_io_particles( name, prt_global_index )
3392
3393	IMPLICIT NONE
3394
3395	CHARACTER(LEN=*), INTENT(IN) :: name !<
3396	INTEGER(idp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: prt_global_index !<
3397
3398	INTEGER(iwp) :: array_size !<
3399	INTEGER(iwp) :: byte_column !<
3400	INTEGER(iwp) :: i !<
3401	INTEGER(iwp) :: ind !<
3402	INTEGER(iwp) :: j !<
3403	INTEGER(iwp) :: k !<
3404	INTEGER(iwp) :: n !<
3405	INTEGER(iwp) :: particle_size !<
3406
3407	INTEGER(KIND=rd_offset_kind) :: disp !<
3408	INTEGER(KIND=rd_offset_kind) :: offset !<
3409	INTEGER(KIND=rd_offset_kind) :: prt_nr_bytes !<
3410
3411	REAL(dp) :: rs !< use REAL*8 to compute max offset
3412	REAL(dp) :: rr !< there is no data type INTEGER*8 in MPI
3413
3414
3415	TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET :: prt_data !<
3416
3417	#if defined( __parallel )
3418	INTEGER, DIMENSION(rd_status_size) :: status !<
3419	#else
3420	TYPE(C_PTR) :: buf
3421	#endif
3422
3423
3424	offset = 0
3425
3426	array_names(header_array_index) = TRIM(name)
3427	array_offset(header_array_index) = array_position
3428	header_array_index = header_array_index + 1
3429
3430	particle_size = STORAGE_SIZE( zero_particle ) / 8
3431
3432	array_size = 0
3433	DO i = nxl, nxr
3434	DO j = nys, nyn
3435	array_size = MAX( array_size, SUM(prt_count(:,j,i)) )
3436	ENDDO
3437	ENDDO
3438
3439	ALLOCATE( prt_data(MAX(array_size,1)) )
3440
3441	!
3442	!-- Write columns of particles.
3443	!-- The particles of a column are stored sequentially in the first dimension of the particles array.
3444	!-- Store only the particles of one cell would be possible with this setup, but the I/O portions
3445	!-- for a maximum of 100 particles are not big enough.
3446	#if defined( __parallel )
3447	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
3448	#endif
3449	prt_nr_bytes = 0
3450	DO i = nxl, nxr
3451	DO j = nys, nyn
3452	disp = array_position + prt_global_index(nzb,j,i) * particle_size
3453	byte_column = SUM( prt_count(:,j,i) ) * particle_size
3454	prt_nr_bytes = MAX( disp+byte_column, prt_nr_bytes )
3455
3456	ind = 1
3457	DO k = nzb, nzt+1
3458	DO n = 1, prt_count(k,j,i)
3459	prt_data(ind) = grid_particles(k,j,i)%particles(n)
3460	ind = ind+1
3461	ENDDO
3462	ENDDO
3463
3464	#if defined( __parallel )
3465	CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
3466	IF ( byte_column > 0 ) THEN
3467	CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
3468	CALL MPI_FILE_WRITE( fh, prt_data, byte_column, MPI_BYTE, status, ierr )
3469	ENDIF
3470	CALL sm_io%sm_node_barrier()
3471	#else
3472	buf = C_LOC(prt_data) ! use C_PTR to avoid another overlay in posix interface
3473	CALL posix_lseek( fh, disp )
3474	CALL posix_write( fh, buf, byte_column )
3475	#endif
3476	ENDDO
3477	ENDDO
3478
3479	#if defined( __parallel )
3480	rs = prt_nr_bytes
3481	CALL MPI_ALLREDUCE( rs, rr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm2d, ierr )
3482	prt_nr_bytes = rr
3483	#else
3484	rr = rs
3485	#endif
3486	array_position = prt_nr_bytes
3487
3488	DEALLOCATE( prt_data )
3489
3490	END SUBROUTINE wrd_mpi_io_particles
3491
3492
3493
3494	!--------------------------------------------------------------------------------------------------!
3495	! Description:
3496	! ------------
3497	!> Write 1d-REAL surface data array with MPI-IO.
3498	!--------------------------------------------------------------------------------------------------!
3499	SUBROUTINE wrd_mpi_io_surface( name, data, first_index )
3500
3501	IMPLICIT NONE
3502
3503	CHARACTER(LEN=*), INTENT(IN) :: name !<
3504
3505	INTEGER(iwp), OPTIONAL :: first_index !<
3506	INTEGER(iwp) :: i !<
3507	INTEGER(iwp) :: j !<
3508	INTEGER(iwp) :: lo_first_index !<
3509	#if defined( __parallel )
3510	INTEGER(iwp) :: buf_start !<
3511	INTEGER(iwp) :: ie !<
3512	INTEGER(iwp) :: is !<
3513	INTEGER(iwp) :: ind_gb !<
3514	INTEGER(iwp) :: ind_out !<
3515	INTEGER(iwp) :: n !<
3516	INTEGER(iwp) :: n_trans !<
3517	#endif
3518
3519	#if defined( __parallel )
3520	INTEGER(KIND=MPI_ADDRESS_KIND) :: disp !< displacement in RMA window
3521	INTEGER(KIND=rd_offset_kind) :: offset !<
3522
3523	INTEGER(iwp), DIMENSION(0:numprocs-1) :: lo_index !<
3524	INTEGER(iwp), DIMENSION(rd_status_size) :: status !<
3525	#endif
3526
3527	REAL(wp), INTENT(IN), DIMENSION(:), TARGET :: data !<
3528	#if defined( __parallel )
3529	REAL(wp), DIMENSION(:), ALLOCATABLE :: get_buffer !<
3530	#endif
3531
3532
3533	lo_first_index = 1
3534
3535	IF ( PRESENT( first_index ) ) THEN
3536	lo_first_index = first_index
3537	ENDIF
3538
3539	!
3540	!-- In case of 2d-data, name is written only once
3541	IF ( lo_first_index == 1 ) THEN
3542
3543	IF ( header_array_index == max_nr_arrays ) THEN
3544	message_string = 'maximum number of 2d/3d-array entries in restart file header exceeded'
3545	CALL message( 'wrd_mpi_io_surface', 'PA0585', 1, 2, 0, 6, 0 )
3546	ENDIF
3547
3548	array_names(header_array_index) = name
3549	array_offset(header_array_index) = array_position
3550	header_array_index = header_array_index + 1
3551
3552	ENDIF
3553
3554	#if defined( __parallel )
3555	offset = 0
3556
3557	ALLOCATE( get_buffer(SUM( transfer_index(4,:) )) )
3558	!
3559	!-- Copy from input array (data) to RMA window to allow the target PEs to get the appropiate data.
3560	!-- At this point, a dummy surface element is added. This makes sure that every x-y gridbox owns
3561	!-- at least one surface element. This way, bookkeeping becomes much easier.
3562	lo_index = nr_surfaces_in_tb
3563	DO i = nxl, nxr
3564	DO j = nys, nyn
3565	is = lo_index(s_index_in_window(j,i)) + 1
3566	ie = is + m_end_index(j,i) - m_start_index(j,i)
3567	!
3568	!-- Store number of surface elements in dummy additional surface element
3569	array_1d(is-1) = e_end_index(j,i) - e_start_index(j,i) + 1
3570	array_1d(is:ie) = data(m_start_index(j,i):m_end_index(j,i))
3571	lo_index(s_index_in_window(j,i)) = lo_index(s_index_in_window(j,i)) + &
3572	e_end_index(j,i) - e_start_index(j,i) + 1
3573	ENDDO
3574	ENDDO
3575	!
3576	!-- On target PE, get data from source PEs which are assigned for output on this PE.
3577	CALL MPI_WIN_FENCE( 0, win_surf, ierr )
3578
3579	buf_start = 1
3580	DO n = 0, numprocs-1
3581	n_trans = transfer_index(4,n)
3582	IF ( n_trans > 0 ) THEN
3583	disp = transfer_index(3,n) - 1
3584	CALL MPI_GET( get_buffer(buf_start), n_trans, MPI_REAL, n, disp, n_trans, MPI_REAL, &
3585	win_surf, ierr )
3586	buf_start = buf_start + n_trans
3587	ENDIF
3588	ENDDO
3589
3590	CALL MPI_WIN_FENCE( 0, win_surf, ierr )
3591	!
3592	!-- Copy data to output buffer. Here, the outpuf buffer matches the indices global_start and
3593	!-- global_end.
3594	ind_gb = 1
3595	DO i = 1, SIZE( local_indices, 2 )
3596	ind_out = local_indices(1,i)
3597	DO j = 1, local_indices(2,i)
3598	array_out(ind_out) = get_buffer(ind_gb)
3599	ind_out = ind_out+1
3600	ind_gb = ind_gb+1
3601	ENDDO
3602	ENDDO
3603
3604	DEALLOCATE( get_buffer )
3605
3606	!
3607	!-- Write data to disk.
3608	CALL sm_io%sm_node_barrier() ! has no effect if I/O on limited number of cores is inactive
3609	IF ( sm_io%iam_io_pe ) THEN
3610	CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_surf, 'native', MPI_INFO_NULL, &
3611	ierr )
3612	CALL MPI_FILE_WRITE_ALL( fh, array_out, SIZE( array_out ), MPI_REAL, status, ierr )
3613	ENDIF
3614	CALL sm_io%sm_node_barrier()
3615	#else
3616	DO i = nxl, nxr
3617	DO j = nys, nyn
3618	array_out(e_start_index(j,i)) = e_end_index(j,i) - e_start_index(j,i) + 1
3619	array_out(e_start_index(j,i)+1:e_end_index(j,i)) = &
3620	data(m_start_index(j,i):m_end_index(j,i))
3621	ENDDO
3622	ENDDO
3623
3624	CALL posix_lseek( fh, array_position )
3625	CALL posix_write( fh, array_out, SIZE(array_out) )
3626	#endif
3627
3628	array_position = array_position + total_number_of_surface_elements * wp
3629
3630	END SUBROUTINE wrd_mpi_io_surface
3631
3632
3633	!--------------------------------------------------------------------------------------------------!
3634	! Description:
3635	! ------------
3636	!> Read 2d-REAL surface data array with MPI-IO.
3637	!> This consist of multiple 1d-REAL surface data arrays.
3638	!--------------------------------------------------------------------------------------------------!
3639	SUBROUTINE wrd_mpi_io_surface_2d( name, data )
3640
3641	IMPLICIT NONE
3642
3643	CHARACTER(LEN=*), INTENT(IN) :: name !<
3644
3645	INTEGER(iwp) :: i !<
3646
3647	REAL(wp), INTENT(IN), DIMENSION(:,:) :: data !<
3648
3649
3650	DO i = 1, SIZE( data, 1 )
3651	CALL wrd_mpi_io_surface( name, data(i,:), first_index = i )
3652	ENDDO
3653
3654	END SUBROUTINE wrd_mpi_io_surface_2d
3655
3656
3657
3658
3659	!--------------------------------------------------------------------------------------------------!
3660	! Description:
3661	! ------------
3662	!> Close restart file for MPI-IO
3663	!--------------------------------------------------------------------------------------------------!
3664	SUBROUTINE rd_mpi_io_close
3665
3666	IMPLICIT NONE
3667
3668	INTEGER(iwp) :: gh_size !<
3669	INTEGER(KIND=rd_offset_kind) :: offset !<
3670	#if defined( __parallel )
3671	INTEGER, DIMENSION(rd_status_size) :: status !<
3672	#endif
3673
3674	#if ! defined( __parallel )
3675	TYPE(C_PTR) :: buf_ptr !<
3676	#endif
3677
3678
3679	offset = 0
3680
3681	IF ( wr_flag .AND. sm_io%iam_io_pe ) THEN
3682
3683	tgh%nr_int = header_int_index - 1
3684	tgh%nr_char = header_char_index - 1
3685	tgh%nr_real = header_real_index - 1
3686	tgh%nr_arrays = header_array_index - 1
3687	tgh%total_nx = iog%nx + 1
3688	tgh%total_ny = iog%ny + 1
3689	tgh%pes_along_x = pdims(1)
3690	tgh%pes_along_y = pdims(2)
3691	IF ( include_total_domain_boundaries ) THEN ! Not sure, if LOGICAL interpretation is the same for all compilers,
3692	tgh%i_outer_bound = 1 ! therefore store as INTEGER in general header
3693	ELSE
3694	tgh%i_outer_bound = 0
3695	ENDIF
3696	!
3697	!-- Check for big/little endian format. This check is currently not used, and could be removed
3698	!-- if we can assume little endian as the default on all machines.
3699	CALL rd_mpi_io_check_endian( tgh%endian )
3700
3701	!
3702	!-- Set default view
3703	#if defined( __parallel )
3704	CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
3705	#endif
3706	!
3707	!-- Write header file
3708	gh_size = storage_size(tgh) / 8
3709	IF ( myid == 0 ) THEN ! myid = 0 always performs I/O, even if I/O is limited to some cores
3710	#if defined( __parallel )
3711	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3712	CALL MPI_FILE_WRITE( fh, tgh, gh_size, MPI_INT, status, ierr )
3713	header_position = header_position + gh_size
3714	!
3715	!-- INTEGER values
3716	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3717	CALL MPI_FILE_WRITE( fh, int_names, SIZE( int_names ) * 32, MPI_CHAR, status, ierr )
3718	header_position = header_position + SIZE( int_names ) * 32
3719
3720	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3721	CALL MPI_FILE_WRITE( fh, int_values, SIZE( int_values ), MPI_INT, status, ierr )
3722	header_position = header_position + SIZE( int_values ) * iwp
3723	!
3724	!-- Character entries
3725	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3726	CALL MPI_FILE_WRITE( fh, text_lines, SIZE( text_lines ) * 128, MPI_CHAR, status, ierr )
3727	header_position = header_position + SIZE( text_lines ) * 128
3728	!
3729	!--- REAL values
3730	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3731	CALL MPI_FILE_WRITE( fh, real_names, SIZE( real_names ) * 32, MPI_CHAR, status, ierr )
3732	header_position = header_position + SIZE( real_names ) * 32
3733
3734	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3735	CALL MPI_FILE_WRITE( fh, real_values, SIZE( real_values ), MPI_REAL, status, ierr )
3736	header_position = header_position + SIZE( real_values ) * wp
3737	!
3738	!-- 2d- and 3d- distributed array headers, all replicated array headers
3739	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3740	CALL MPI_FILE_WRITE( fh, array_names, SIZE( array_names ) * 32, MPI_CHAR, status, ierr )
3741	header_position = header_position + SIZE( array_names ) * 32
3742
3743	CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
3744	CALL MPI_FILE_WRITE( fh, array_offset, SIZE( array_offset ) * MPI_OFFSET_KIND, MPI_BYTE, &
3745	status, ierr ) ! There is no I*8 datatype in Fortran
3746	header_position = header_position + SIZE( array_offset ) * rd_offset_kind
3747	#else
3748	CALL posix_lseek( fh, header_position )
3749	buf_ptr = C_LOC( tgh )
3750	CALL posix_write( fh, buf_ptr, gh_size )
3751	header_position = header_position + gh_size
3752	!
3753	!-- INTEGER values
3754	CALL posix_lseek( fh, header_position )
3755	CALL posix_write( fh, int_names )
3756	header_position = header_position + SIZE( int_names ) * 32
3757
3758	CALL posix_lseek( fh, header_position )
3759	CALL posix_write( fh, int_values, SIZE( int_values ) )
3760	header_position = header_position + SIZE( int_values ) * iwp
3761	!
3762	!-- Character entries
3763	CALL posix_lseek( fh, header_position )
3764	CALL posix_write( fh, text_lines )
3765	header_position = header_position + SIZE( text_lines ) * 128
3766	!
3767	!-- REAL values
3768	CALL posix_lseek( fh, header_position )
3769	CALL posix_write( fh, real_names )
3770	header_position = header_position + SIZE( real_names ) * 32
3771
3772	CALL posix_lseek( fh, header_position )
3773	CALL posix_write( fh, real_values, SIZE( real_values ) )
3774	header_position = header_position + SIZE( real_values ) * wp
3775	!
3776	!-- 2d- and 3d-distributed array headers, all replicated array headers
3777	CALL posix_lseek( fh, header_position )
3778	CALL posix_write( fh, array_names )
3779	header_position = header_position + SIZE( array_names ) * 32
3780
3781	CALL posix_lseek( fh, header_position )
3782	CALL posix_write( fh, array_offset, SIZE( array_offset ) )
3783	header_position = header_position + SIZE( array_offset ) * rd_offset_kind
3784	#endif
3785	IF ( debug_output ) CALL rd_mpi_io_print_header
3786	ENDIF
3787
3788	ENDIF
3789
3790	!
3791	!-- Free file types
3792	CALL rd_mpi_io_free_filetypes
3793
3794	!
3795	!-- Close MPI-IO files
3796	#if defined( __parallel )
3797	!
3798	!-- Free RMA windows
3799	IF ( cyclic_fill_mode ) THEN
3800	CALL MPI_WIN_FREE( rmawin_2di, ierr )
3801	CALL MPI_WIN_FREE( rmawin_2d, ierr )
3802	CALL MPI_WIN_FREE( rmawin_3d, ierr )
3803	ENDIF
3804	#endif
3805
3806	IF ( ALLOCATED( e_start_index ) ) DEALLOCATE( e_start_index )
3807	IF ( ALLOCATED( e_end_index ) ) DEALLOCATE( e_end_index )
3808	IF ( ALLOCATED( m_start_index ) ) DEALLOCATE( m_start_index )
3809	IF ( ALLOCATED( m_end_index ) ) DEALLOCATE( m_end_index )
3810	IF ( ALLOCATED( nr_surfaces_in_tb ) ) DEALLOCATE( nr_surfaces_in_tb )
3811	IF ( ALLOCATED( s_index_in_tb ) ) DEALLOCATE( s_index_in_tb )
3812	IF ( ALLOCATED( s_index_in_window ) ) DEALLOCATE( s_index_in_window )
3813	IF ( ALLOCATED( transfer_index ) ) DEALLOCATE( transfer_index )
3814
3815	IF ( .NOT. pe_active_for_read ) RETURN
3816	!
3817	!-- TODO: better explain the following message
3818	!-- In case on non cyclic read, pe_active_for_read is set .TRUE.
3819	IF ( sm_io%iam_io_pe ) THEN
3820
3821	#if defined( __parallel )
3822	CALL MPI_FILE_CLOSE( fh, ierr )
3823	#else
3824	CALL posix_close( fh )
3825	#endif
3826
3827	ENDIF
3828
3829	fh = -1
3830
3831	restart_file_size = array_position / ( 1024.0_dp * 1024.0_dp )
3832
3833	END SUBROUTINE rd_mpi_io_close
3834
3835
3836
3837	FUNCTION rd_mpi_io_check_open() RESULT ( isopen )
3838
3839	LOGICAL :: isopen
3840
3841	isopen = ( fh /= -1 )
3842
3843	END FUNCTION rd_mpi_io_check_open
3844
3845
3846
3847	!--------------------------------------------------------------------------------------------------!
3848	! Description:
3849	! ------------
3850	!> This subroutine prepares a filetype and some variables for the I/O of surface data.
3851	!> Whenever a new set of start_index and end_index is used, rd_mpi_io_surface_filetypes has to be
3852	!> called. A main feature of this subroutine is computing the global start indices of the 1d- and
3853	!> 2d- surface arrays.
3854	!> Even if I/O is done by a limited number of cores only, the surface data are read by ALL cores!
3855	!> Reading them by some cores and then distributing the data would result in complicated code
3856	!> which is suspectable for errors and overloads the reading subroutine. Since reading of surface
3857	!> data is not time critical (data size is comparably small), it will be read by all cores.
3858	!--------------------------------------------------------------------------------------------------!
3859	RECURSIVE SUBROUTINE rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, &
3860	global_start, global_end )
3861
3862	IMPLICIT NONE
3863
3864	INTEGER(iwp) :: e_lo_start !<
3865	INTEGER(iwp) :: i !< loop index
3866	INTEGER(iwp) :: j !< loop index
3867	INTEGER(iwp) :: index_offset !<
3868	INTEGER(iwp) :: last_end_index !<
3869	INTEGER(iwp) :: lo_start !<
3870	INTEGER(iwp) :: nr_surfcells_pe !<
3871	INTEGER(iwp) :: rest_cells_pe !<
3872	INTEGER(iwp) :: rest_bound !<
3873	#if defined( __parallel )
3874	INTEGER(iwp) :: io_end_index !<
3875	INTEGER(iwp) :: io_start_index !<
3876	INTEGER(iwp) :: n !< loop index
3877	INTEGER(iwp) :: nr_previous !<
3878	#endif
3879
3880	INTEGER(iwp), DIMENSION(0:numprocs-1,2) :: nr_surfcells_all_s !<
3881	INTEGER(iwp), DIMENSION(0:numprocs-1,2) :: nr_surfcells_all_r !<
3882	#if defined( __parallel )
3883	INTEGER(iwp), DIMENSION(1) :: dims1 !< global dimension for MPI_TYPE_CREATE_SUBARRAY
3884	INTEGER(iwp), DIMENSION(1) :: lsize1 !< local size for MPI_TYPE_CREATE_SUBARRAY
3885	INTEGER(iwp), DIMENSION(0:numprocs-1) :: nr_gp_with_surfaces_for_pe !< number of horizontal gridpoints containing surface elements for the respective pe
3886	INTEGER(iwp), DIMENSION(0:numprocs-1) :: nr_surface_elements_for_pe !< total number of surface elements for the respective pe
3887	INTEGER(iwp), DIMENSION(0:pdims(1)) :: nr_surf_cells_x !<
3888	INTEGER(iwp), DIMENSION(0:pdims(1)) :: nr_surf_cells_x_s !<
3889	INTEGER(iwp), DIMENSION(1) :: start1 !< start index for MPI_TYPE_CREATE_SUBARRAY
3890	INTEGER(iwp), DIMENSION(nxl:nxr) :: sum_y !<
3891	#endif
3892
3893	INTEGER(iwp), INTENT(INOUT), DIMENSION(nys:nyn,nxl:nxr) :: end_index !< local end indx
3894	INTEGER(iwp), INTENT(INOUT), DIMENSION(nys:nyn,nxl:nxr) :: global_start !< global start index
3895	INTEGER(iwp), INTENT(INOUT), DIMENSION(nys:nyn,nxl:nxr) :: global_end !< global end index
3896	INTEGER(iwp), INTENT(INOUT), DIMENSION(nys:nyn,nxl:nxr) :: start_index !< local start index
3897	#if defined( __parallel )
3898	INTEGER(iwp), DIMENSION(0:myidy,nxl:nxr) :: nr_previous_y !<
3899	INTEGER(iwp), DIMENSION(0:pdims(2),nxl:nxr) :: nr_surf_cells_y !<
3900	INTEGER(iwp), DIMENSION(0:pdims(2),nxl:nxr) :: nr_surf_cells_y_s !<
3901	INTEGER(iwp), DIMENSION(4,0:numprocs-1) :: transfer_index_s !<
3902	#endif
3903
3904	LOGICAL, INTENT(OUT) :: data_to_write !< returns .TRUE., if surface data have been written
3905	LOGICAL :: only_dummy_elements !< only dummy elements, i.e. no data to write
3906
3907
3908	!
3909	!-- Clean up previous calls.
3910	#if defined( __parallel )
3911	IF ( win_surf /= -1 ) THEN
3912	CALL MPI_WIN_FREE( win_surf, ierr )
3913	DEALLOCATE( array_1d )
3914	win_surf = -1
3915	ENDIF
3916	IF ( ft_surf /= -1 .AND. sm_io%iam_io_pe ) THEN
3917	CALL MPI_TYPE_FREE( ft_surf, ierr )
3918	ENDIF
3919	ft_surf = -1
3920	IF ( sm_io%is_sm_active() ) THEN
3921	IF ( win_out /= -1 ) THEN
3922	CALL MPI_WIN_FREE( win_out, ierr )
3923	win_out = -1
3924	ENDIF
3925	ELSE
3926	IF ( ASSOCIATED( array_out ) ) DEALLOCATE( array_out )
3927	ENDIF
3928	#else
3929	IF ( ASSOCIATED( array_out ) ) DEALLOCATE( array_out )
3930	#endif
3931
3932	IF ( cyclic_fill_mode ) THEN
3933	CALL cyclic_fill_surface_filetype
3934	RETURN
3935	ELSE
3936	IF ( .NOT. ALLOCATED( e_end_index ) ) ALLOCATE( e_end_index(nys:nyn,nxl:nxr) )
3937	IF ( .NOT. ALLOCATED( e_start_index ) ) ALLOCATE( e_start_index(nys:nyn,nxl:nxr) )
3938	IF ( .NOT. ALLOCATED( m_end_index ) ) ALLOCATE( m_end_index(nys:nyn,nxl:nxr) )
3939	IF ( .NOT. ALLOCATED( m_start_index ) ) ALLOCATE( m_start_index(nys:nyn,nxl:nxr) )
3940	IF ( .NOT. ALLOCATED( nr_surfaces_in_tb ) ) ALLOCATE( nr_surfaces_in_tb(0:numprocs-1) )
3941	IF ( .NOT. ALLOCATED( s_index_in_tb ) ) ALLOCATE( s_index_in_tb(0:numprocs-1) )
3942	IF ( .NOT. ALLOCATED( s_index_in_window ) ) ALLOCATE( s_index_in_window(nys:nyn,nxl:nxr) )
3943	IF ( .NOT. ALLOCATED( transfer_index ) ) ALLOCATE( transfer_index(4,0:numprocs-1) )
3944	ENDIF
3945
3946	IF ( wr_flag) THEN
3947	!
3948	!-- Add one dummy value at every grid box.
3949	!-- This allows to use MPI_FILE_WRITE_ALL and MPI_FILE_READ_ALL with subarray file type.
3950	index_offset = 0
3951	last_end_index = 0
3952	DO i = nxl, nxr
3953	DO j = nys, nyn
3954	e_start_index(j,i) = start_index (j,i) + index_offset
3955	IF ( end_index (j,i) - start_index(j,i) < 0 ) THEN
3956	e_end_index (j,i) = last_end_index+1
3957	last_end_index = last_end_index+1
3958	ELSE
3959	e_end_index (j,i) = end_index(j,i) + index_offset + 1
3960	last_end_index = e_end_index (j,i)
3961	ENDIF
3962	index_offset = index_offset + 1
3963	ENDDO
3964	ENDDO
3965	#if defined( __parallel )
3966	!
3967	!-- Compute indices for global, PE independent 1-d surface element array.
3968	nr_surf_cells_y_s = 0
3969	!
3970	!-- Count number of surface elements in y-direction for every x.
3971	DO i = nxl, nxr
3972	nr_surf_cells_y_s(myidy,i) = SUM( e_end_index (:,i) - e_start_index (:,i) + 1 )
3973	ENDDO
3974	!
3975	!-- Distribute these elements to all PEs along y.
3976	CALL MPI_ALLREDUCE( nr_surf_cells_y_s, nr_surf_cells_y, SIZE( nr_surf_cells_y ), &
3977	MPI_INTEGER, MPI_SUM, comm1dy, ierr )
3978	!
3979	!-- Sum all surface elements along y for individual x PEs
3980	nr_surf_cells_x_s = 0
3981	nr_surf_cells_x_s(myidx) = SUM( nr_surf_cells_y )
3982	!
3983	!-- Distribute to all PEs along x.
3984	CALL MPI_ALLREDUCE( nr_surf_cells_x_s, nr_surf_cells_x, SIZE( nr_surf_cells_x ), &
3985	MPI_INTEGER, MPI_SUM, comm1dx, ierr )
3986	DO i = nxl, nxr
3987	nr_previous_y(:,i) = 0
3988	DO n = 1, myidy
3989	nr_previous_y(n,i) = nr_previous_y(n-1,i) + nr_surf_cells_y(n-1,i)
3990	ENDDO
3991	ENDDO
3992
3993	sum_y(nxl) = SUM( nr_surf_cells_y(:,nxl) )
3994	DO i = nxl, nxr
3995	IF ( i > nxl ) THEN
3996	sum_y(i) = sum_y(i-1) + SUM( nr_surf_cells_y(:,i) )
3997	ENDIF
3998	ENDDO
3999
4000	nr_previous = 0
4001	IF ( myidx >= 1 ) THEN
4002	nr_previous = SUM(nr_surf_cells_x(0:myidx-1))
4003	ENDIF
4004
4005	global_start(nys,nxl) = 1 + nr_previous + nr_previous_y(myidy,nxl)
4006	DO j = nys+1, nyn
4007	global_start(j,nxl) = global_start(j-1,nxl) + e_end_index(j-1,nxl) - &
4008	e_start_index(j-1,nxl) + 1
4009	ENDDO
4010
4011	DO i = nxl+1, nxr
4012	global_start(nys,i) = 1 + nr_previous + nr_previous_y(myidy,i) + sum_y(i-1)
4013	DO j = nys+1, nyn
4014	global_start(j,i) = global_start(j-1,i) + e_end_index(j-1,i) - e_start_index(j-1,i) + 1
4015	ENDDO
4016	ENDDO
4017	#else
4018	global_start = e_start_index
4019	#endif
4020	DO i = nxl, nxr
4021	DO j = nys, nyn
4022	global_end(j,i) = global_start(j,i) + e_end_index (j,i) - e_start_index (j,i)
4023	ENDDO
4024	ENDDO
4025
4026	ELSE
4027	!
4028	!-- In case of read, compute e_start_index and e_end_index for current processor grid.
4029	!-- This data contains one extra value for every i and j.
4030	e_lo_start = 1
4031	lo_start = 1
4032	DO i = nxl, nxr
4033	DO j = nys, nyn
4034	e_start_index(j,i) = e_lo_start
4035	e_end_index(j,i) = e_lo_start + global_end(j,i) - global_start(j,i)
4036	e_lo_start = e_lo_start + global_end(j,i) - global_start(j,i) + 1
4037	start_index(j,i) = lo_start
4038	end_index(j,i) = lo_start + global_end(j,i) - global_start(j,i) - 1
4039	lo_start = lo_start + global_end(j,i) - global_start(j,i)
4040	ENDDO
4041	ENDDO
4042
4043	ENDIF
4044
4045	nr_surfcells_all_s = 0
4046	nr_surfcells_all_s(myid,1) = MAXVAL( e_end_index ) ! don't split surface elements of one gridbox
4047	nr_surfcells_all_s(myid,2) = MAXVAL( e_end_index - e_start_index )
4048
4049	#if defined( __parallel )
4050	CALL MPI_ALLREDUCE( nr_surfcells_all_s, nr_surfcells_all_r, SIZE( nr_surfcells_all_s ), &
4051	MPI_INTEGER, MPI_SUM, comm2d, ierr )
4052	#else
4053	nr_surfcells_all_r = nr_surfcells_all_s
4054	#endif
4055
4056	total_number_of_surface_elements = 0
4057	DO i = 0, numprocs-1
4058	IF ( i == myid ) THEN
4059	glo_start = total_number_of_surface_elements + 1
4060	ENDIF
4061	total_number_of_surface_elements = total_number_of_surface_elements + nr_surfcells_all_r(i,1)
4062	ENDDO
4063	only_dummy_elements = ( MAXVAL( nr_surfcells_all_r(:,2) ) <= 0 )
4064
4065	!
4066	!-- Compute indices of equally distributed surface elements.
4067	!-- Number of surface elements scheduled for ouput on this PE:
4068	nr_surfcells_pe = total_number_of_surface_elements / numprocs
4069	rest_cells_pe = MOD( total_number_of_surface_elements, numprocs )
4070	rest_bound = rest_cells_pe * ( nr_surfcells_pe + 1 )
4071	m_start_index = start_index
4072	m_end_index = end_index
4073
4074	!
4075	!-- Compute number of surface elements on source PE, which have to be send to the corresponding
4076	!-- target PE.
4077	#if defined( __parallel )
4078	nr_gp_with_surfaces_for_pe = 0
4079	nr_surface_elements_for_pe = 0
4080	DO i = nxl, nxr
4081	DO j = nys, nyn
4082	IF ( rest_cells_pe == 0 ) THEN
4083	s_index_in_window(j,i) = ( global_start(j,i) - 1 ) / nr_surfcells_pe
4084	ELSE
4085	IF ( global_start(j,i) <= rest_bound ) THEN
4086	s_index_in_window(j,i) = ( global_start(j,i) - 1 ) / ( nr_surfcells_pe + 1 )
4087	ELSE
4088	s_index_in_window(j,i) = ( global_start(j,i) - rest_bound - 1 ) / nr_surfcells_pe
4089	s_index_in_window(j,i) = s_index_in_window(j,i) + rest_cells_pe
4090	ENDIF
4091	!
4092	!-- TODO: Test output, to be removed later.
4093	IF ( s_index_in_window(j,i) >= numprocs ) THEN
4094	WRITE( 9,'(A,8I8)' ) 's_index_in_window ', j, i, s_index_in_window(j,i), &
4095	global_start(j,i) , nr_surfcells_pe
4096	FLUSH( 9 )
4097	CALL MPI_ABORT( comm2d, 1, ierr )
4098	ENDIF
4099	ENDIF
4100	nr_gp_with_surfaces_for_pe(s_index_in_window(j,i)) = &
4101	nr_gp_with_surfaces_for_pe(s_index_in_window(j,i)) + 1
4102	nr_surface_elements_for_pe(s_index_in_window(j,i)) = &
4103	nr_surface_elements_for_pe(s_index_in_window(j,i)) + &
4104	e_end_index(j,i) - e_start_index(j,i) + 1
4105	ENDDO
4106	ENDDO
4107
4108	!
4109	!-- Compute start index in the transfer buffer on the source side for the corresponding target PE.
4110	s_index_in_tb(0) = 1
4111	nr_surfaces_in_tb(0) = 1
4112	DO n = 1, numprocs-1
4113	s_index_in_tb(n) = s_index_in_tb(n-1) + nr_gp_with_surfaces_for_pe(n-1)
4114	nr_surfaces_in_tb(n) = nr_surfaces_in_tb(n-1) + nr_surface_elements_for_pe(n-1)
4115	ENDDO
4116	!
4117	!-- Buffer distribution on the source side.
4118	DO n = 0, numprocs-1
4119	transfer_index_s(1,n) = s_index_in_tb(n)
4120	transfer_index_s(2,n) = nr_gp_with_surfaces_for_pe(n)
4121	transfer_index_s(3,n) = nr_surfaces_in_tb(n)
4122	transfer_index_s(4,n) = nr_surface_elements_for_pe(n)
4123	ENDDO
4124
4125	CALL MPI_ALLTOALL( transfer_index_s, 4, MPI_INTEGER, transfer_index, 4, MPI_INTEGER, comm2d, &
4126	ierr)
4127	!
4128	!-- Buffer distribution on the target side side.
4129	CALL get_remote_indices()
4130	!
4131	!-- Create surface element file type.
4132	IF ( total_number_of_surface_elements > 0 .AND. .NOT. only_dummy_elements) THEN
4133	data_to_write = .TRUE.
4134	ELSE
4135	data_to_write = .FALSE.
4136	ENDIF
4137
4138	CALL MPI_ALLREDUCE( global_end(nyn,nxr), dims1(1), 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
4139	start1(1) = MINVAL( local_indices(1,:) ) - 1
4140	IF ( sm_io%is_sm_active() ) THEN
4141	CALL MPI_ALLREDUCE( SUM( local_indices(2,:) ), lsize1(1), 1, MPI_INTEGER, MPI_SUM, &
4142	sm_io%comm_shared, ierr )
4143	ELSE
4144	lsize1(1) = SUM( local_indices(2,:) )
4145	ENDIF
4146
4147	IF ( sm_io%iam_io_pe ) THEN
4148	IF ( total_number_of_surface_elements > 0 ) THEN
4149	CALL MPI_TYPE_CREATE_SUBARRAY( 1, dims1, lsize1, start1, MPI_ORDER_FORTRAN, MPI_REAL, &
4150	ft_surf, ierr )
4151	CALL MPI_TYPE_COMMIT( ft_surf, ierr )
4152	ENDIF
4153	ENDIF
4154	!
4155	!-- Allocate rma window to supply surface data to other PEs.
4156	CALL rd_alloc_rma_mem( array_1d, SUM( nr_surface_elements_for_pe ), win_surf )
4157	!
4158	!-- Allocate shared array on IO-PE to supply data for MPI-IO (write or read).
4159	IF ( sm_io%is_sm_active() ) THEN
4160	IF ( sm_io%iam_io_pe ) THEN
4161	io_start_index = start1(1) + 1
4162	io_end_index = start1(1) + lsize1(1)
4163	ENDIF
4164	CALL MPI_BCAST( io_start_index, 1, MPI_INTEGER, 0, sm_io%comm_shared, ierr )
4165	CALL MPI_BCAST( io_end_index, 1, MPI_INTEGER, 0, sm_io%comm_shared, ierr )
4166	CALL sm_io%sm_allocate_shared( array_out, io_start_index, io_end_index, win_out )
4167	ELSE
4168	ALLOCATE( array_out(start1(1)+1:start1(1)+lsize1(1)) )
4169	ENDIF
4170	#else
4171	IF ( total_number_of_surface_elements > 0 .AND. .NOT. only_dummy_elements ) THEN
4172	data_to_write = .TRUE.
4173	ELSE
4174	data_to_write = .FALSE.
4175	ENDIF
4176	ALLOCATE( array_out(1:total_number_of_surface_elements) )
4177	#endif
4178
4179	CONTAINS
4180
4181	SUBROUTINE cyclic_fill_surface_filetype
4182
4183	INTEGER(iwp) :: i !< loop index
4184	INTEGER(iwp) :: j !< loop index
4185
4186
4187	IF ( .NOT. ALLOCATED( o_start_index ) ) ALLOCATE( o_start_index(nys:nyn,nxl:nxr) )
4188	IF ( .NOT. ALLOCATED( o_end_index ) ) ALLOCATE( o_end_index(nys:nyn,nxl:nxr) )
4189
4190	lo_start = 1
4191	DO i = nxl, nxr
4192	DO j = nys, nyn
4193	o_start_index(j,i) = lo_start
4194	o_end_index(j,i) = lo_start
4195	lo_start = lo_start + 1
4196	ENDDO
4197	ENDDO
4198	start_index = o_start_index
4199	end_index = o_end_index
4200
4201	IF ( MAXVAL( global_end-global_start ) > 1 ) THEN
4202	message_string = 'cylic-fill method does not allow more than one surface element ' // &
4203	'per grid box'
4204	CALL message( 'cyclic_fill_surface_filetype', 'PA0742', 1, 2, 0, 6, 0 )
4205	ENDIF
4206	!
4207	!-- Activate grid of the smaller prerun, i.e. grid variables like nxl, nxr, nys, nyn and others
4208	!-- are set according to the prerun layout.
4209	CALL prerun_grid%activate_grid_from_this_class()
4210
4211	IF ( pe_active_for_read ) THEN
4212
4213	IF ( .NOT. ALLOCATED( c_global_start ) ) ALLOCATE( c_global_start(nys:nyn,nxl:nxr) )
4214	IF ( .NOT. ALLOCATED( c_global_end ) ) ALLOCATE( c_global_end(nys:nyn,nxl:nxr) )
4215	IF ( .NOT. ALLOCATED( c_start_index ) ) ALLOCATE( c_start_index(nys:nyn,nxl:nxr) )
4216	IF ( .NOT. ALLOCATED( c_end_index ) ) ALLOCATE( c_end_index(nys:nyn,nxl:nxr) )
4217
4218	DO i = nxl, nxr
4219	DO j = nys, nyn
4220	c_global_start(j,i) = global_start(j,i)
4221	c_global_end(j,i) = global_end(j,i)
4222	ENDDO
4223	ENDDO
4224	!
4225	!-- Recursive call of rd_mpi_io_surface_filetypes.
4226	!-- Prerun data are read, but they are treated as if they are mainrun data, just on a smaller
4227	!-- grid.
4228	cyclic_fill_mode = .FALSE.
4229	CALL rd_mpi_io_surface_filetypes( c_start_index, c_end_index, data_to_write, &
4230	c_global_start, c_global_end )
4231	cyclic_fill_mode = .TRUE.
4232
4233	ENDIF
4234	!
4235	!-- Activate grid of the mainrun, i.e. grid variables like nxl, nxr, nys, nyn and others
4236	!-- are set according to the mainrun layout.
4237	CALL mainrun_grid%activate_grid_from_this_class()
4238
4239	#if defined( __parallel )
4240	CALL MPI_BCAST( data_to_write, 1, MPI_LOGICAL, 0, comm2d, ierr )
4241	#endif
4242
4243	END SUBROUTINE cyclic_fill_surface_filetype
4244
4245	#if defined( __parallel )
4246	!
4247	!-- Get the indices of the surface elements inside the RMA window on the remote PE.
4248	!-- This information is required to fetch the surface element data on remote PEs
4249	!-- in rrd_mpi_io_surface and wrd_mpi_io_surface.
4250	SUBROUTINE get_remote_indices
4251
4252	INTEGER(iwp) :: buf_start !<
4253	INTEGER(iwp) :: i !<
4254	INTEGER(iwp) :: j !<
4255	INTEGER(iwp) :: n !<
4256	INTEGER(iwp) :: n_trans !<
4257	INTEGER(iwp) :: win_ind !<
4258
4259	INTEGER(KIND=MPI_ADDRESS_KIND) :: disp !< displacement in RMA window
4260	INTEGER(KIND=MPI_ADDRESS_KIND) :: winsize !< size of RMA window
4261
4262	INTEGER(iwp), DIMENSION(0:numprocs-1) :: lo_index !<
4263
4264	INTEGER(iwp), POINTER, DIMENSION(:,:) :: surf_val_index !<
4265
4266
4267	IF ( ALLOCATED( local_indices ) ) DEALLOCATE( local_indices )
4268	ALLOCATE( local_indices(2,MAX( SUM( transfer_index(2,:) ), 2 )))
4269
4270	local_indices(1,:) = 0
4271	local_indices(2,:) = 0
4272
4273	winsize = MAX( 2 * SUM( nr_gp_with_surfaces_for_pe ), 2 )
4274
4275	ALLOCATE( surf_val_index(2,winsize) )
4276	winsize = winsize * iwp
4277	CALL MPI_WIN_CREATE( surf_val_index, winsize, iwp, MPI_INFO_NULL, comm2d, win_ind, ierr )
4278	CALL MPI_WIN_FENCE( 0, win_ind, ierr )
4279
4280	lo_index = s_index_in_tb
4281	DO i = nxl, nxr
4282	DO j = nys, nyn
4283	surf_val_index(1,lo_index(s_index_in_window(j,i))) = global_start(j,i)
4284	surf_val_index(2,lo_index(s_index_in_window(j,i))) = global_end(j,i) - &
4285	global_start(j,i) + 1
4286	lo_index(s_index_in_window(j,i)) = lo_index(s_index_in_window(j,i)) + 1
4287	ENDDO
4288	ENDDO
4289
4290	CALL MPI_WIN_FENCE( 0, win_ind, ierr )
4291
4292	buf_start = 1
4293	DO n = 0, numprocs-1
4294	n_trans = transfer_index(2,n)
4295	IF ( n_trans > 0 ) THEN
4296	disp = 2 * ( transfer_index(1,n) - 1 )
4297	CALL MPI_GET( local_indices(1,buf_start), 2n_trans, MPI_INTEGER, n, disp, 2n_trans, &
4298	MPI_INTEGER, win_ind, ierr )
4299	buf_start = buf_start + n_trans
4300	ENDIF
4301	ENDDO
4302
4303	CALL MPI_WIN_FENCE( 0, win_ind, ierr )
4304
4305	buf_start = 1
4306	DO n = 0, numprocs-1
4307	n_trans = transfer_index(2,n)
4308	IF ( n_trans > 0 ) THEN
4309	disp = transfer_index(1,n) - 1
4310	buf_start = buf_start + n_trans
4311	ENDIF
4312	ENDDO
4313
4314	CALL MPI_WIN_FREE( win_ind, ierr )
4315
4316	DEALLOCATE( surf_val_index )
4317
4318	END SUBROUTINE get_remote_indices
4319
4320	!--------------------------------------------------------------------------------------------------!
4321	! Description:
4322	! ------------
4323	!> Allocate memory and create window for one-sided communication (1-d INTEGER array)
4324	!--------------------------------------------------------------------------------------------------!
4325	SUBROUTINE rd_alloc_rma_mem( array, idim, win )
4326
4327	IMPLICIT NONE
4328
4329	INTEGER(iwp), INTENT(IN) :: idim !< Dimension of this 1-D array
4330	INTEGER :: ierr !< MPI error code
4331	INTEGER(iwp), INTENT(OUT) :: win !< MPI window
4332	INTEGER(KIND=MPI_ADDRESS_KIND) :: winsize !< size of RMA window
4333
4334	REAL(wp), DIMENSION(:), POINTER, INTENT(INOUT) :: array !< array to access RMA window locally
4335
4336
4337	winsize = MAX( idim, 2 )
4338	ALLOCATE( array(winsize) )
4339	winsize = winsize * wp
4340	CALL MPI_WIN_CREATE( array, winsize, wp, MPI_INFO_NULL, comm2d, win, ierr )
4341	array = -1
4342	CALL MPI_WIN_FENCE( 0, win, ierr )
4343
4344	END SUBROUTINE rd_alloc_rma_mem
4345	#endif
4346
4347	END SUBROUTINE rd_mpi_io_surface_filetypes
4348
4349
4350	!--------------------------------------------------------------------------------------------------!
4351	! Description:
4352	! ------------
4353	!> This subroutine creates file types to access 2d-/3d-REAL arrays and 2d-INTEGER arrays
4354	!> distributed in blocks among processes to a single file that contains the global arrays.
4355	!--------------------------------------------------------------------------------------------------!
4356	SUBROUTINE rd_mpi_io_create_filetypes
4357
4358	IMPLICIT NONE
4359
4360	INTEGER, DIMENSION(2) :: dims2 !<
4361	INTEGER, DIMENSION(2) :: lize2 !<
4362	INTEGER, DIMENSION(2) :: start2 !<
4363
4364	INTEGER, DIMENSION(3) :: dims3 !<
4365	INTEGER, DIMENSION(3) :: lize3 !<
4366	INTEGER, DIMENSION(3) :: start3 !<
4367
4368	TYPE(domain_decomposition_grid_features) :: save_io_grid !< temporary variable to store grid settings
4369
4370
4371	IF ( sm_io%is_sm_active() ) THEN
4372	save_io_grid = sm_io%io_grid
4373	ENDIF
4374
4375	IF( .NOT. pe_active_for_read ) RETURN
4376
4377	IF ( cyclic_fill_mode ) CALL prerun_grid%activate_grid_from_this_class()
4378
4379	!
4380	!-- Decision, if storage with or without ghost layers.
4381	!-- Please note that the indexing of the global array always starts at 0, even in Fortran.
4382	!-- Therefore the PE local indices have to be shifted by nbgp in the case with ghost layers.
4383	IF ( include_total_domain_boundaries ) THEN
4384
4385	iog%nxl = nxl + nbgp
4386	iog%nxr = nxr + nbgp
4387	iog%nys = nys + nbgp
4388	iog%nyn = nyn + nbgp
4389	iog%nnx = nnx
4390	iog%nny = nny
4391	iog%nx = nx + 2 * nbgp
4392	iog%ny = ny + 2 * nbgp
4393	IF ( myidx == 0 ) THEN
4394	iog%nxl = iog%nxl - nbgp
4395	iog%nnx = iog%nnx + nbgp
4396	ENDIF
4397	IF ( myidx == pdims(1)-1 ) THEN
4398	iog%nxr = iog%nxr + nbgp
4399	iog%nnx = iog%nnx + nbgp
4400	ENDIF
4401	IF ( myidy == 0 ) THEN
4402	iog%nys = iog%nys - nbgp
4403	iog%nny = iog%nny + nbgp
4404	ENDIF
4405	IF ( myidy == pdims(2)-1 ) THEN
4406	iog%nyn = iog%nyn + nbgp
4407	iog%nny = iog%nny + nbgp
4408	ENDIF
4409
4410	CALL sm_io%sm_adjust_outer_boundary()
4411
4412	ELSE
4413
4414	iog%nxl = nxl
4415	iog%nxr = nxr
4416	iog%nys = nys
4417	iog%nyn = nyn
4418	iog%nnx = nnx
4419	iog%nny = nny
4420	iog%nx = nx
4421	iog%ny = ny
4422
4423	ENDIF
4424
4425	IF ( sm_io%is_sm_active() ) THEN
4426	#if defined( __parallel )
4427	CALL sm_io%sm_allocate_shared( array_2d, sm_io%io_grid%nxl, sm_io%io_grid%nxr, &
4428	sm_io%io_grid%nys, sm_io%io_grid%nyn, win_2dr )
4429	CALL sm_io%sm_allocate_shared( array_2di, save_io_grid%nxl, save_io_grid%nxr, &
4430	save_io_grid%nys, save_io_grid%nyn, win_2di )
4431	CALL sm_io%sm_allocate_shared( array_3d, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr, &
4432	sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3dr )
4433	CALL sm_io%sm_allocate_shared( array_3di4, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,&
4434	sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3di4 )
4435	#endif
4436	ELSE
4437	ALLOCATE( array_2d(iog%nxl:iog%nxr,iog%nys:iog%nyn) )
4438	ALLOCATE( array_2di(nxl:nxr,nys:nyn) )
4439	ALLOCATE( array_3d(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
4440	ALLOCATE( array_3di4(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
4441	sm_io%io_grid = iog
4442	ENDIF
4443
4444	!
4445	!-- Create filetype for 2d-REAL array with ghost layers around the total domain
4446	dims2(1) = iog%nx + 1
4447	dims2(2) = iog%ny + 1
4448
4449	lize2(1) = sm_io%io_grid%nnx
4450	lize2(2) = sm_io%io_grid%nny
4451
4452	start2(1) = sm_io%io_grid%nxl
4453	start2(2) = sm_io%io_grid%nys
4454
4455	#if defined( __parallel )
4456	IF ( sm_io%iam_io_pe ) THEN
4457	CALL MPI_TYPE_CREATE_SUBARRAY( 2, dims2, lize2, start2, MPI_ORDER_FORTRAN, MPI_REAL, &
4458	ft_2d, ierr )
4459	CALL MPI_TYPE_COMMIT( ft_2d, ierr )
4460	ENDIF
4461	#endif
4462	!
4463	!-- Create filetype for 2d-INTEGER array without ghost layers around the total domain
4464	dims2(1) = nx + 1
4465	dims2(2) = ny + 1
4466
4467	IF ( sm_io%is_sm_active() ) THEN
4468
4469	lize2(1) = save_io_grid%nnx
4470	lize2(2) = save_io_grid%nny
4471
4472	start2(1) = save_io_grid%nxl
4473	start2(2) = save_io_grid%nys
4474
4475	ELSE
4476
4477	lize2(1) = nnx
4478	lize2(2) = nny
4479
4480	start2(1) = nxl
4481	start2(2) = nys
4482
4483	ENDIF
4484
4485	#if defined( __parallel )
4486	IF ( sm_io%iam_io_pe ) THEN
4487	CALL MPI_TYPE_CREATE_SUBARRAY( 2, dims2, lize2, start2, MPI_ORDER_FORTRAN, MPI_INTEGER, &
4488	ft_2di_nb, ierr )
4489	CALL MPI_TYPE_COMMIT( ft_2di_nb, ierr )
4490	ENDIF
4491	#endif
4492	!
4493	!-- Create filetype for 3d-REAL array
4494	dims3(1) = nz + 2
4495	dims3(2) = iog%nx + 1
4496	dims3(3) = iog%ny + 1
4497
4498	lize3(1) = dims3(1)
4499	lize3(2) = sm_io%io_grid%nnx
4500	lize3(3) = sm_io%io_grid%nny
4501
4502	start3(1) = nzb
4503	start3(2) = sm_io%io_grid%nxl
4504	start3(3) = sm_io%io_grid%nys
4505
4506	#if defined( __parallel )
4507	IF ( sm_io%iam_io_pe ) THEN
4508	CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_REAL, ft_3d, &
4509	ierr )
4510	CALL MPI_TYPE_COMMIT( ft_3d, ierr )
4511	ENDIF
4512	#endif
4513
4514	IF ( cyclic_fill_mode ) CALL mainrun_grid%activate_grid_from_this_class()
4515
4516	END SUBROUTINE rd_mpi_io_create_filetypes
4517
4518
4519
4520	!--------------------------------------------------------------------------------------------------!
4521	! Description:
4522	! ------------
4523	!> This subroutine creates file types to access 3d-INTEGER4 arrays and 3d-INTEGER8 arrays
4524	!> distributed in blocks among processes to a single file that contains the global arrays.
4525	!> These types are only used for particle data.
4526	!--------------------------------------------------------------------------------------------------!
4527	SUBROUTINE rd_mpi_io_particle_filetypes
4528
4529	IMPLICIT NONE
4530
4531	INTEGER, DIMENSION(3) :: dims3 !<
4532	INTEGER, DIMENSION(3) :: lize3 !<
4533	INTEGER, DIMENSION(3) :: start3 !<
4534
4535	TYPE(domain_decomposition_grid_features) :: save_io_grid !< temporary variable to store grid settings
4536
4537	!
4538	!-- MPI_INTEGER8 is not standard MPI, but is supported on most MPI distibutions
4539	!-- If not suppported, a workaround could be enabled with the following preprocessor directive
4540	!#if defined( __NO_INTEGER8)
4541	! CALL MPI_TYPE_CONTIGUOUS( 2, MPI_INTEGER, MPI_INTEGER8, ierr )
4542	! CALL MPI_TYPE_COMMIT( MPI_INTEGER8, ierr )
4543	!#endif
4544
4545	IF ( sm_io%is_sm_active() ) THEN
4546	save_io_grid = sm_io%io_grid
4547	ENDIF
4548
4549	IF( .NOT. pe_active_for_read ) RETURN
4550
4551	!
4552	!-- Decision, if storage with or without ghost layers.
4553	!-- Please note that the indexing of the global array always starts at 0, even in Fortran.
4554	!-- Therefore the PE local indices have to be shifted by nbgp in the case with ghost layers.
4555	IF ( include_total_domain_boundaries ) THEN
4556
4557	iog%nxl = nxl + nbgp
4558	iog%nxr = nxr + nbgp
4559	iog%nys = nys + nbgp
4560	iog%nyn = nyn + nbgp
4561	iog%nnx = nnx
4562	iog%nny = nny
4563	iog%nx = nx + 2 * nbgp
4564	iog%ny = ny + 2 * nbgp
4565	IF ( myidx == 0 ) THEN
4566	iog%nxl = iog%nxl - nbgp
4567	iog%nnx = iog%nnx + nbgp
4568	ENDIF
4569	IF ( myidx == pdims(1)-1 ) THEN
4570	iog%nxr = iog%nxr + nbgp
4571	iog%nnx = iog%nnx + nbgp
4572	ENDIF
4573	IF ( myidy == 0 ) THEN
4574	iog%nys = iog%nys - nbgp
4575	iog%nny = iog%nny + nbgp
4576	ENDIF
4577	IF ( myidy == pdims(2)-1 ) THEN
4578	iog%nyn = iog%nyn + nbgp
4579	iog%nny = iog%nny + nbgp
4580	ENDIF
4581
4582	CALL sm_io%sm_adjust_outer_boundary()
4583
4584	ELSE
4585
4586	iog%nxl = nxl
4587	iog%nxr = nxr
4588	iog%nys = nys
4589	iog%nyn = nyn
4590	iog%nnx = nnx
4591	iog%nny = nny
4592	iog%nx = nx
4593	iog%ny = ny
4594
4595	ENDIF
4596
4597	IF ( sm_io%is_sm_active() ) THEN
4598	#if defined( __parallel )
4599	CALL sm_io%sm_allocate_shared( array_3di8, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,&
4600	sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3di8 )
4601	#endif
4602	ELSE
4603	ALLOCATE( array_3di8(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
4604
4605	sm_io%io_grid = iog
4606	ENDIF
4607
4608	!
4609	!-- Create filetype for 3d INTEGER array
4610	dims3(1) = nz + 2
4611	dims3(2) = iog%nx + 1
4612	dims3(3) = iog%ny + 1
4613
4614	lize3(1) = dims3(1)
4615	lize3(2) = sm_io%io_grid%nnx
4616	lize3(3) = sm_io%io_grid%nny
4617
4618	start3(1) = nzb
4619	start3(2) = sm_io%io_grid%nxl
4620	start3(3) = sm_io%io_grid%nys
4621
4622	#if defined( __parallel )
4623	IF ( sm_io%iam_io_pe ) THEN
4624	CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_INTEGER, &
4625	ft_3di4, ierr )
4626	CALL MPI_TYPE_COMMIT( ft_3di4, ierr )
4627
4628	CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_INTEGER8, &
4629	ft_3di8, ierr )
4630	CALL MPI_TYPE_COMMIT( ft_3di8, ierr )
4631	ENDIF
4632	#endif
4633
4634	END SUBROUTINE rd_mpi_io_particle_filetypes
4635
4636
4637
4638	!--------------------------------------------------------------------------------------------------!
4639	! Description:
4640	! ------------
4641	!> This subroutine creates file types to access 3d-soil arrays distributed in blocks among processes
4642	!> to a single file that contains the global arrays. It is not required for the serial mode.
4643	!--------------------------------------------------------------------------------------------------!
4644	SUBROUTINE rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
4645
4646	IMPLICIT NONE
4647
4648	INTEGER, INTENT(IN) :: nzb_soil !<
4649	INTEGER, INTENT(IN) :: nzt_soil !<
4650
4651	#if defined( __parallel )
4652	INTEGER, DIMENSION(3) :: dims3 !<
4653	INTEGER, DIMENSION(3) :: lize3 !<
4654	INTEGER, DIMENSION(3) :: start3 !<
4655
4656
4657	IF ( sm_io%is_sm_active() ) THEN
4658	CALL sm_io%sm_allocate_shared( array_3d_soil, nzb_soil, nzt_soil, sm_io%io_grid%nxl, &
4659	sm_io%io_grid%nxr, sm_io%io_grid%nys, sm_io%io_grid%nyn, &
4660	win_3ds )
4661	ELSE
4662	ALLOCATE( array_3d_soil(nzb_soil:nzt_soil,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
4663	sm_io%io_grid = iog
4664	ENDIF
4665
4666	!
4667	!-- Create filetype for 3d-soil array
4668	dims3(1) = nzt_soil - nzb_soil + 1
4669	dims3(2) = iog%nx + 1
4670	dims3(3) = iog%ny + 1
4671
4672	lize3(1) = dims3(1)
4673	lize3(2) = sm_io%io_grid%nnx
4674	lize3(3) = sm_io%io_grid%nny
4675
4676	start3(1) = nzb_soil
4677	start3(2) = sm_io%io_grid%nxl
4678	start3(3) = sm_io%io_grid%nys
4679
4680	IF ( sm_io%iam_io_pe ) THEN
4681	CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_REAL, &
4682	ft_3dsoil, ierr )
4683	CALL MPI_TYPE_COMMIT( ft_3dsoil, ierr )
4684	ENDIF
4685	#else
4686	ALLOCATE( array_3d_soil(nzb_soil:nzt_soil,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
4687	sm_io%io_grid = iog
4688	#endif
4689
4690	END SUBROUTINE rd_mpi_io_create_filetypes_3dsoil
4691
4692	!--------------------------------------------------------------------------------------------------!
4693	! Description:
4694	! ------------
4695	!> Free all file types that have been created for MPI-IO.
4696	!--------------------------------------------------------------------------------------------------!
4697	SUBROUTINE rd_mpi_io_free_filetypes
4698
4699	IMPLICIT NONE
4700
4701	#if defined( __parallel )
4702	IF ( filetypes_created ) THEN
4703
4704	IF ( sm_io%iam_io_pe ) THEN
4705	CALL MPI_TYPE_FREE( ft_2d, ierr )
4706	CALL MPI_TYPE_FREE( ft_2di_nb, ierr )
4707	CALL MPI_TYPE_FREE( ft_3d, ierr )
4708	ENDIF
4709
4710	IF ( sm_io%is_sm_active() ) THEN
4711	CALL sm_io%sm_free_shared( win_2dr )
4712	CALL sm_io%sm_free_shared( win_2di )
4713	CALL sm_io%sm_free_shared( win_3dr )
4714	ELSE
4715	DEALLOCATE( array_2d, array_2di, array_3d )
4716	ENDIF
4717
4718	ENDIF
4719
4720	!
4721	!-- Free last surface filetype
4722	IF ( sm_io%iam_io_pe .AND. ft_surf /= -1 ) THEN
4723	CALL MPI_TYPE_FREE( ft_surf, ierr )
4724	ENDIF
4725
4726	IF ( sm_io%is_sm_active() .AND. win_surf /= -1 ) THEN
4727	CALL sm_io%sm_free_shared( win_surf )
4728	ENDIF
4729
4730	!
4731	!-- Free last particle filetypes
4732	IF ( sm_io%iam_io_pe .AND. ft_3di4 /= -1 ) THEN
4733	CALL MPI_TYPE_FREE( ft_3di4, ierr )
4734	ft_3di4 = -1
4735	ENDIF
4736	IF ( sm_io%iam_io_pe .AND. ft_3di8 /= -1 ) THEN
4737	CALL MPI_TYPE_FREE( ft_3di8, ierr )
4738	ft_3di8 = -1
4739	ENDIF
4740
4741	IF ( sm_io%is_sm_active() .AND. win_3di4 /= -1 ) THEN
4742	CALL sm_io%sm_free_shared( win_3di4 )
4743	win_3di4 = -1
4744	ENDIF
4745	IF ( sm_io%is_sm_active() .AND. win_3di8 /= -1 ) THEN
4746	CALL sm_io%sm_free_shared( win_3di8 )
4747	win_3di8 = -1
4748	ENDIF
4749
4750	IF ( win_start /= -1 ) THEN
4751	CALL sm_io%sm_free_shared( win_start)
4752	CALL sm_io%sm_free_shared( win_end)
4753	CALL sm_io%sm_free_shared( win_glost)
4754	win_start = -1
4755	win_end = -1
4756	win_glost = -1
4757	ENDIF
4758
4759	ft_surf = -1
4760	win_surf = -1
4761	#else
4762	IF ( ASSOCIATED( array_2d ) ) DEALLOCATE( array_2d )
4763	IF ( ASSOCIATED( array_2di ) ) DEALLOCATE( array_2di )
4764	IF ( ASSOCIATED( array_3d ) ) DEALLOCATE( array_3d )
4765	IF ( ASSOCIATED( array_3di4 ) ) DEALLOCATE( array_3di4 )
4766	IF ( ASSOCIATED( array_3di8 ) ) DEALLOCATE( array_3di8 )
4767	#endif
4768
4769	END SUBROUTINE rd_mpi_io_free_filetypes
4770
4771
4772
4773	!--------------------------------------------------------------------------------------------------!
4774	! Description:
4775	! ------------
4776	!> Print the restart data file header (MPI-IO format) for debugging.
4777	!--------------------------------------------------------------------------------------------------!
4778	SUBROUTINE rd_mpi_io_print_header
4779
4780	IMPLICIT NONE
4781
4782	INTEGER(iwp) :: i !<
4783
4784
4785	WRITE (9,*) 'header position after reading the restart file header: ', header_position
4786	WRITE (9,*) ' '
4787	WRITE (9,*) 'restart file header content:'
4788	WRITE (9,*) '----------------------------'
4789	WRITE (9,*) ' '
4790
4791	WRITE (9,*) ' CHARACTER header values Total number: ', tgh%nr_char
4792	WRITE (9,*) ' '
4793	DO i = 1, tgh%nr_char
4794	WRITE( 9, '(I3,A,1X,A)' ) i, ': ', text_lines(i)(1:80)
4795	ENDDO
4796	WRITE (9,*) ' '
4797
4798	WRITE (9,*) ' INTEGER header variables and values Total number: ', tgh%nr_int
4799	WRITE (9,*) ' '
4800	DO i = 1, tgh%nr_int
4801	WRITE(9,*) ' variable: ', int_names(i), ' value: ', int_values(i)
4802	ENDDO
4803	WRITE (9,*) ' '
4804
4805	WRITE (9,*) ' REAL header variables and values Total number: ', tgh%nr_real
4806	WRITE (9,*) ' '
4807	DO i = 1, tgh%nr_real
4808	WRITE(9,*) ' variable: ', real_names(i), ' value: ', real_values(i)
4809	ENDDO
4810	WRITE (9,*) ' '
4811
4812	WRITE (9,*) ' Header entries with offset (2d/3d arrays) Total number: ', tgh%nr_arrays
4813	WRITE (9,*) ' '
4814	DO i = 1, tgh%nr_arrays
4815	WRITE(9,*) ' variable: ', array_names(i), ' offset: ', array_offset(i)
4816	ENDDO
4817	WRITE (9,*) ' '
4818
4819	END SUBROUTINE rd_mpi_io_print_header
4820
4821
4822
4823	!--------------------------------------------------------------------------------------------------!
4824	! Description:
4825	! ------------
4826	!> Check if big/little endian data format is used.
4827	!> An int4 pointer is set to a int8 variable, the int*8 is set to 1, and then it is checked, if
4828	!> the first 4 bytes of the pointer are equal 1 (little endian) or not.
4829	!--------------------------------------------------------------------------------------------------!
4830	SUBROUTINE rd_mpi_io_check_endian( i_endian )
4831
4832	IMPLICIT NONE
4833
4834	INTEGER, INTENT(out) :: i_endian !<
4835	INTEGER(KIND=8), TARGET :: int8 !<
4836
4837	INTEGER, DIMENSION(1) :: bufshape !<
4838	INTEGER(KIND=4), POINTER, DIMENSION(:) :: int4 !<
4839
4840	TYPE(C_PTR) :: ptr !<
4841
4842
4843	ptr = C_LOC( int8 )
4844	bufshape(1) = 2
4845	CALL C_F_POINTER( ptr, int4, bufshape )
4846
4847	int8 = 1
4848
4849	IF ( int4(1) == 1 ) THEN
4850	i_endian = 1 ! Little endian
4851	ELSE
4852	i_endian = 2 ! Big endian
4853	ENDIF
4854
4855	END SUBROUTINE rd_mpi_io_check_endian
4856
4857	END MODULE restart_data_mpi_io_mod

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

Context Navigation

source: palm/trunk/SOURCE/restart_data_mpi_io_mod.f90

Download in other formats: