Home

Context Navigation

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

Last change on this file since 4830 was 4828, checked in by Giersch, 4 years ago
Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code
Property svn:keywords set to `Id`
File size: 152.0 KB

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