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

Last change on this file since 4811 was 4735, checked in by raasch, 4 years ago

interface for 3d logical arrays added

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