Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

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

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