Home

Context Navigation

source: palm/trunk/SOURCE/restart_data_mpi_io_mod.f90 @ 4893

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |