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

Last change on this file since 4621 was 4621, checked in by raasch, 4 months ago

further bugfix for r4618: mising preprocessor directives for serial mode added

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