source: palm/trunk/SOURCE/restart_data_mpi_io_mod.f90

Last change on this file was 4896, checked in by raasch, 2 months ago

small re-formatting to follow the coding standard, typo in file appendix removed, more meaningful variable names assigned, redundant code removed

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