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

Last change on this file since 4631 was 4629, checked in by raasch, 4 years ago

support for MPI Fortran77 interface (mpif.h) removed

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