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

Last change on this file since 4712 was 4694, checked in by pavelkrc, 4 years ago

Fix writing and reading of surface data to/from MPI restart file + small fixes

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