Home

Context Navigation

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

Last change on this file since 4716 was 4694, checked in by pavelkrc, 4 years ago
Fix writing and reading of surface data to/from MPI restart file + small fixes
Property svn:keywords set to `Id`
File size: 150.1 KB

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