Home

Context Navigation

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

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