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

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

bugfix for r4618: unused variable removed

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