Home

Context Navigation

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

Last change on this file since 4834 was 4832, checked in by raasch, 4 years ago
bugfix in redblack algorithm: lower i,j indices need to start alternatively with even or odd value on the coarsest grid level, if the subdomain has an uneven number of gridpoints along x/y; some debug output flushed
Property svn:keywords set to `Id`
File size: 152.1 KB

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