Home

Context Navigation

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

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