Home

Context Navigation

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

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