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

Last change on this file since 4591 was 4591, checked in by raasch, 8 months ago

files re-formatted to follow the PALM coding standard

  • Property svn:keywords set to Id
File size: 99.5 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 4591 2020-07-06 15:56:08Z raasch $
27! File re-formatted to follow the PALM coding standard
28!
29!
30! 4539 2020-05-18 14:05:17Z raasch
31! Checks added, if index limits in header are exceeded
32! Bugfix in rrd_mpi_io_int_2d
33!
34! 4536 2020-05-17 17:24:13Z raasch
35! Messages and debug output converted to PALM routines
36!
37! 4534 2020-05-14 18:35:22Z raasch
38! I/O on reduced number of cores added (using shared memory MPI)
39!
40! 4500 2020-04-17 10:12:45Z suehring
41! Fix too long lines
42!
43! 4498 2020-04-15 14:26:31Z raasch
44! Bugfix for creation of filetypes, argument removed from rd_mpi_io_open
45!
46! 4497 2020-04-15 10:20:51Z raasch
47! Last bugfix deactivated because of compile problems
48!
49! 4496 2020-04-15 08:37:26Z raasch
50! Problem with posix read arguments for surface data fixed
51!
52! 4495 2020-04-13 20:11:20Z raasch
53! Initial version (K. Ketelsen), adjusted to PALM formatting standards (s. Raasch)
54!
55!
56!
57! Description:
58! ------------
59!> Routines for restart data handling using MPI-IO.
60!--------------------------------------------------------------------------------------------------!
61 MODULE restart_data_mpi_io_mod
62
63#if defined( __parallel )
64#if defined( __mpifh )
65    INCLUDE "mpif.h"
66#else
67    USE MPI
68#endif
69#else
70    USE posix_interface,                                                                           &
71        ONLY:  posix_close,                                                                        &
72               posix_lseek,                                                                        &
73               posix_open,                                                                         &
74               posix_read,                                                                         &
75               posix_write
76#endif
77
78    USE, INTRINSIC ::  ISO_C_BINDING
79
80    USE control_parameters,                                                                        &
81        ONLY:  debug_output,                                                                       &
82               debug_string,                                                                       &
83               include_total_domain_boundaries,                                                    &
84               message_string,                                                                     &
85               restart_data_format_input,                                                          &
86               restart_data_format_output,                                                         &
87               restart_file_size
88
89    USE exchange_horiz_mod,                                                                        &
90        ONLY:  exchange_horiz,                                                                     &
91               exchange_horiz_2d
92
93    USE indices,                                                                                   &
94        ONLY:  nbgp,                                                                               &
95               nnx,                                                                                &
96               nny,                                                                                &
97               nx,                                                                                 &
98               nxl,                                                                                &
99               nxlg,                                                                               &
100               nxr,                                                                                &
101               nxrg,                                                                               &
102               ny,                                                                                 &
103               nyn,                                                                                &
104               nyng,                                                                               &
105               nys,                                                                                &
106               nysg,                                                                               &
107               nz,                                                                                 &
108               nzb,                                                                                &
109               nzt
110
111    USE kinds
112
113    USE pegrid,                                                                                    &
114        ONLY:  comm1dx,                                                                            &
115               comm1dy,                                                                            &
116               comm2d,                                                                             &
117               myid,                                                                               &
118               myidx,                                                                              &
119               myidy,                                                                              &
120               npex,                                                                               &
121               npey,                                                                               &
122               numprocs,                                                                           &
123               pdims
124
125    USE shared_memory_io_mod,                                                                      &
126        ONLY:  local_boundaries,                                                                   &
127               sm_class
128
129
130    IMPLICIT NONE
131
132    CHARACTER(LEN=128) ::  io_file_name  !> internal variable to communicate filename between different subroutines
133
134#if defined( __parallel )
135    INTEGER(iwp)            ::  ierr                              !< error status of MPI-calls
136    INTEGER(iwp), PARAMETER ::  rd_offset_kind = MPI_OFFSET_KIND  !< Adress or Offset kind
137    INTEGER(iwp), PARAMETER ::  rd_status_size = MPI_STATUS_SIZE  !<
138#else
139    INTEGER(iwp), PARAMETER ::  rd_offset_kind = C_SIZE_T         !<
140    INTEGER(iwp), PARAMETER ::  rd_status_size = 1                !< Not required in sequential mode
141#endif
142
143    INTEGER(iwp)            ::  debug_level = 1  !< TODO: replace with standard debug output steering
144
145    INTEGER(iwp)            ::  comm_io          !< Communicator for MPI-IO
146    INTEGER(iwp)            ::  fh               !< MPI-IO file handle
147#if defined( __parallel )
148    INTEGER(iwp)            ::  fhs = -1         !< MPI-IO file handle to open file with comm2d always
149#endif
150    INTEGER(iwp)            ::  ft_surf = -1     !< MPI filetype surface data
151#if defined( __parallel )
152    INTEGER(iwp)            ::  ft_2di_nb        !< MPI filetype 2D array INTEGER no outer boundary
153    INTEGER(iwp)            ::  ft_2d            !< MPI filetype 2D array REAL with outer boundaries
154    INTEGER(iwp)            ::  ft_3d            !< MPI filetype 3D array REAL with outer boundaries
155    INTEGER(iwp)            ::  ft_3dsoil        !< MPI filetype for 3d-soil array
156#endif
157    INTEGER(iwp)            ::  glo_start        !< global start index on this PE
158#if defined( __parallel )
159    INTEGER(iwp)            ::  local_start      !<
160#endif
161    INTEGER(iwp)            ::  nr_iope          !<
162    INTEGER(iwp)            ::  nr_val           !< local number of values in x and y direction
163#if defined( __parallel )
164    INTEGER(iwp)            ::  win_2di          !<
165    INTEGER(iwp)            ::  win_2dr          !<
166    INTEGER(iwp)            ::  win_3dr          !<
167    INTEGER(iwp)            ::  win_3ds          !<
168    INTEGER(iwp)            ::  win_surf = -1    !<
169#endif
170    INTEGER(iwp)            ::  total_number_of_surface_values  !< total number of values for one variable
171
172    INTEGER(KIND=rd_offset_kind) ::  array_position   !<
173    INTEGER(KIND=rd_offset_kind) ::  header_position  !<
174
175    INTEGER(iwp), DIMENSION(:,:), POINTER, CONTIGUOUS ::  array_2di  !<
176
177    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_end_index     !<
178    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_global_start  !<
179    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_start_index   !<
180
181
182    LOGICAL ::  all_pes_write                 !< all PEs have data to write
183    LOGICAL ::  filetypes_created             !<
184    LOGICAL ::  io_on_limited_cores_per_node  !< switch to shared memory MPI-IO
185    LOGICAL ::  rd_flag                       !< file is opened for read
186    LOGICAL ::  wr_flag                       !< file is opened for write
187
188#if defined( __parallel )
189    REAL(wp), DIMENSION(:), POINTER, CONTIGUOUS     ::  array_1d       !<
190#endif
191    REAL(wp), DIMENSION(:,:), POINTER, CONTIGUOUS   ::  array_2d       !<
192    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  array_3d       !<
193    REAL(wp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  array_3d_soil  !<
194
195!
196!-- Handling of outer boundaries
197    TYPE(local_boundaries) ::  lb  !<
198
199!
200!-- General Header (first 32 byte in restart file)
201    TYPE general_header
202       INTEGER(iwp) :: endian         !< little endian (1) or big endian (2) internal format
203       INTEGER(iwp) :: i_outer_bound  !< if 1, outer boundaries are stored in restart file
204       INTEGER(iwp) :: nr_arrays      !< number of arrays in restart files
205       INTEGER(iwp) :: nr_char        !< number of Text strings entries in header
206       INTEGER(iwp) :: nr_int         !< number of INTEGER entries in header
207       INTEGER(iwp) :: nr_real        !< number of REAL entries in header
208       INTEGER(iwp) :: total_nx       !< total number of points in x-direction
209       INTEGER(iwp) :: total_ny       !< total number of points in y-direction
210    END TYPE general_header
211
212    TYPE(general_header), TARGET ::  tgh    !<
213
214    TYPE(sm_class)               ::  sm_io  !<
215
216!
217!-- Declaration of varibales for file header section
218    INTEGER(KIND=rd_offset_kind)                ::  header_int_index
219    INTEGER, PARAMETER                          ::  max_nr_int=256
220    CHARACTER(LEN=32), DIMENSION(max_nr_int)    ::  int_names
221    INTEGER(KIND=iwp), DIMENSION(max_nr_int)    ::  int_values
222
223    INTEGER(KIND=rd_offset_kind)                ::  header_char_index
224    INTEGER, PARAMETER                          ::  max_nr_char=128
225    CHARACTER(LEN=128), DIMENSION(max_nr_char)  ::  text_lines
226
227    INTEGER(KIND=rd_offset_kind)                ::  header_real_index
228    INTEGER, PARAMETER                          ::  max_nr_real=256
229    CHARACTER(LEN=32), DIMENSION(max_nr_real)   ::  real_names
230    REAL(KIND=wp), DIMENSION(max_nr_real)       ::  real_values
231
232    INTEGER(KIND=rd_offset_kind)                ::  header_array_index
233    INTEGER, PARAMETER                          ::  max_nr_arrays=600
234    CHARACTER(LEN=32), DIMENSION(max_nr_arrays) ::  array_names
235    INTEGER(KIND=rd_offset_kind), DIMENSION(max_nr_arrays) :: array_offset
236    SAVE
237
238    PRIVATE
239
240    PUBLIC  restart_file_size, total_number_of_surface_values
241
242!
243!-- PALM interfaces
244    INTERFACE rd_mpi_io_check_array
245       MODULE PROCEDURE rd_mpi_io_check_array
246    END INTERFACE rd_mpi_io_check_array
247
248    INTERFACE rd_mpi_io_close
249       MODULE PROCEDURE rd_mpi_io_close
250    END INTERFACE rd_mpi_io_close
251
252    INTERFACE rd_mpi_io_open
253       MODULE PROCEDURE rd_mpi_io_open
254    END INTERFACE rd_mpi_io_open
255
256    INTERFACE rrd_mpi_io
257       MODULE PROCEDURE rrd_mpi_io_char
258       MODULE PROCEDURE rrd_mpi_io_int
259       MODULE PROCEDURE rrd_mpi_io_int_2d
260       MODULE PROCEDURE rrd_mpi_io_logical
261       MODULE PROCEDURE rrd_mpi_io_real
262       MODULE PROCEDURE rrd_mpi_io_real_2d
263       MODULE PROCEDURE rrd_mpi_io_real_3d
264       MODULE PROCEDURE rrd_mpi_io_real_3d_soil
265    END INTERFACE rrd_mpi_io
266
267    INTERFACE rrd_mpi_io_global_array
268       MODULE PROCEDURE rrd_mpi_io_global_array_int_1d
269       MODULE PROCEDURE rrd_mpi_io_global_array_real_1d
270       MODULE PROCEDURE rrd_mpi_io_global_array_real_2d
271       MODULE PROCEDURE rrd_mpi_io_global_array_real_3d
272       MODULE PROCEDURE rrd_mpi_io_global_array_real_4d
273    END INTERFACE rrd_mpi_io_global_array
274
275    INTERFACE rrd_mpi_io_surface
276       MODULE PROCEDURE rrd_mpi_io_surface
277       MODULE PROCEDURE rrd_mpi_io_surface_2d
278    END INTERFACE rrd_mpi_io_surface
279
280    INTERFACE rd_mpi_io_surface_filetypes
281       MODULE PROCEDURE rd_mpi_io_surface_filetypes
282    END INTERFACE rd_mpi_io_surface_filetypes
283
284    INTERFACE wrd_mpi_io
285       MODULE PROCEDURE wrd_mpi_io_char
286       MODULE PROCEDURE wrd_mpi_io_int
287       MODULE PROCEDURE wrd_mpi_io_int_2d
288       MODULE PROCEDURE wrd_mpi_io_logical
289       MODULE PROCEDURE wrd_mpi_io_real
290       MODULE PROCEDURE wrd_mpi_io_real_2d
291       MODULE PROCEDURE wrd_mpi_io_real_3d
292       MODULE PROCEDURE wrd_mpi_io_real_3d_soil
293    END INTERFACE wrd_mpi_io
294
295    INTERFACE wrd_mpi_io_global_array
296       MODULE PROCEDURE wrd_mpi_io_global_array_int_1d
297       MODULE PROCEDURE wrd_mpi_io_global_array_real_1d
298       MODULE PROCEDURE wrd_mpi_io_global_array_real_2d
299       MODULE PROCEDURE wrd_mpi_io_global_array_real_3d
300       MODULE PROCEDURE wrd_mpi_io_global_array_real_4d
301    END INTERFACE wrd_mpi_io_global_array
302
303    INTERFACE wrd_mpi_io_surface
304       MODULE PROCEDURE wrd_mpi_io_surface
305       MODULE PROCEDURE wrd_mpi_io_surface_2d
306    END INTERFACE wrd_mpi_io_surface
307
308    PUBLIC  rd_mpi_io_check_array,                                                                 &
309            rd_mpi_io_close,                                                                       &
310            rd_mpi_io_open,                                                                        &
311            rrd_mpi_io,                                                                            &
312            rrd_mpi_io_global_array,                                                               &
313            rrd_mpi_io_surface,                                                                    &
314            rd_mpi_io_surface_filetypes,                                                           &
315            wrd_mpi_io,                                                                            &
316            wrd_mpi_io_global_array,                                                               &
317            wrd_mpi_io_surface
318
319
320 CONTAINS
321
322
323!--------------------------------------------------------------------------------------------------!
324! Description:
325! ------------
326!> Open restart file for read or write with MPI-IO
327!--------------------------------------------------------------------------------------------------!
328 SUBROUTINE rd_mpi_io_open( action, file_name, open_for_global_io_only )
329
330    IMPLICIT NONE
331
332    CHARACTER(LEN=*), INTENT(IN) ::  action     !<
333    CHARACTER(LEN=*), INTENT(IN) ::  file_name  !<
334
335    INTEGER(iwp)                 ::  i          !<
336    INTEGER(iwp)                 ::  gh_size    !<
337
338    INTEGER(KIND=rd_offset_kind) ::  offset     !<
339
340#if defined( __parallel )
341    INTEGER, DIMENSION(rd_status_size) ::  status  !<
342#endif
343
344    LOGICAL, INTENT(IN), OPTIONAL ::  open_for_global_io_only  !<
345    LOGICAL                       ::  set_filetype             !<
346
347#if ! defined( __parallel )
348    TYPE(C_PTR)                   ::  buf_ptr  !<
349#endif
350
351
352    offset = 0
353    io_on_limited_cores_per_node = .FALSE.
354
355    rd_flag = ( TRIM( action ) == 'READ'  .OR. TRIM( action ) == 'read'  )
356    wr_flag = ( TRIM( action ) == 'WRITE' .OR. TRIM( action ) == 'write' )
357
358    IF ( .NOT. ( rd_flag .OR. wr_flag ) )  THEN
359       message_string = 'illegal action "' // TRIM( action ) // '" for opening restart files'
360       CALL message( 'restart_data_mpi_io_mod', 'PA0720', 1, 2, 0, 6, 0 )
361    ENDIF
362!
363!-- Store name of I/O file to communicate it internally within this module.
364    io_file_name = file_name
365!
366!-- Setup for IO on a limited number of threads per node (using shared memory MPI)
367    IF ( rd_flag )  THEN
368       set_filetype = .TRUE.
369       IF ( TRIM( restart_data_format_input ) == 'mpi_shared_memory' )  THEN
370          io_on_limited_cores_per_node = .TRUE.
371       ENDIF
372    ENDIF
373
374    IF ( TRIM( restart_data_format_output ) == 'mpi_shared_memory' .AND.  wr_flag )  THEN
375       io_on_limited_cores_per_node = .TRUE.
376    ENDIF
377!
378!-- Shared memory MPI is not used for reading of global data
379    IF ( PRESENT( open_for_global_io_only )  .AND.  rd_flag )  THEN
380       IF ( open_for_global_io_only )  THEN
381          io_on_limited_cores_per_node = .FALSE.
382          set_filetype                 = .FALSE.
383       ENDIF
384    ENDIF
385
386    CALL sm_io%sm_init_comm( io_on_limited_cores_per_node )
387
388!
389!-- Set communicator to be used. If all cores are doing I/O, comm2d is used as usual.
390    IF( sm_io%is_sm_active() )  THEN
391       comm_io = sm_io%comm_io
392    ELSE
393       comm_io = comm2d
394    ENDIF
395
396!
397!-- Create subarrays and file types
398    filetypes_created = .FALSE.
399
400!
401!-- In case of read it is not known yet if data include total domain. Filetypes will be created
402!-- further below.
403    IF ( wr_flag )  THEN
404       CALL rd_mpi_io_create_filetypes
405       filetypes_created = .TRUE.
406    ENDIF
407
408!
409!-- Open file for MPI-IO
410#if defined( __parallel )
411    IF ( sm_io%iam_io_pe )  THEN
412
413       IF ( rd_flag )  THEN
414
415          IF ( debug_output )  THEN
416             WRITE( debug_string, * )  'open joint restart file "' // TRIM( io_file_name ) //      &
417                                       '" for read with MPI-IO'
418             CALL debug_message( debug_string, 'start' )
419          ENDIF
420
421          CALL MPI_FILE_OPEN( comm_io, TRIM( io_file_name ), MPI_MODE_RDONLY, MPI_INFO_NULL, fh,   &
422                              ierr )
423
424          IF ( ierr /= 0 )  THEN
425             message_string = 'error opening restart file "' // TRIM( io_file_name ) //            &
426                              '" for reading with MPI-IO'
427             CALL message( 'rrd_mpi_io_open', 'PA0727', 3, 2, 0, 6, 0 )
428          ENDIF
429
430          IF ( debug_output )  THEN
431             WRITE( debug_string, * )  'open joint restart file "' // TRIM( io_file_name ) //      &
432                                       '" for read with MPI-IO'
433             CALL debug_message( debug_string, 'end' )
434          ENDIF
435
436       ELSEIF ( wr_flag )  THEN
437
438          IF ( debug_output )  THEN
439             WRITE( debug_string, * )  'open joint restart file "' // TRIM( io_file_name ) //      &
440                                       '" for write with MPI-IO'
441             CALL debug_message( debug_string, 'start' )
442          ENDIF
443
444          CALL MPI_FILE_OPEN( comm_io, TRIM( io_file_name ), MPI_MODE_CREATE+MPI_MODE_WRONLY,      &
445                              MPI_INFO_NULL, fh, ierr )
446
447          IF ( ierr /= 0 )  THEN
448             message_string = 'error opening restart file "' // TRIM( io_file_name ) //            &
449                              '" for writing with MPI-IO'
450             CALL message( 'rrd_mpi_io_open', 'PA0728', 3, 2, 0, 6, 0 )
451          ENDIF
452
453          IF ( debug_output )  THEN
454             WRITE( debug_string, * )  'open joint restart file "' // TRIM( io_file_name ) //      &
455                                       '" for write with MPI-IO'
456             CALL debug_message( debug_string, 'end' )
457          ENDIF
458
459       ENDIF
460
461    ENDIF
462#else
463    IF ( rd_flag )  THEN
464
465       IF ( debug_output )  THEN
466          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
467                                    '" for read in serial mode (posix)'
468          CALL debug_message( debug_string, 'start' )
469       ENDIF
470
471       fh = posix_open( TRIM( io_file_name ), .TRUE. )
472
473       IF ( debug_output )  THEN
474          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
475                                    '" for read in serial mode (posix)'
476          CALL debug_message( debug_string, 'end' )
477       ENDIF
478
479    ELSEIF ( wr_flag )  THEN
480
481       IF ( debug_output )  THEN
482          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
483                                    '" for write in serial mode (posix)'
484          CALL debug_message( debug_string, 'start' )
485       ENDIF
486
487       fh = posix_open( TRIM( io_file_name ), .FALSE. )
488
489       IF ( debug_output )  THEN
490          WRITE( debug_string, * )  'open restart file "' // TRIM( io_file_name ) //               &
491                                    '" for write in serial mode (posix)'
492          CALL debug_message( debug_string, 'end' )
493       ENDIF
494
495    ENDIF
496
497    IF ( fh < 0 )  THEN
498       message_string = 'error opening restart file for posix I/O'
499       CALL message( 'restart_data_mpi_io_mod', 'PA0721', 1, 2, 0, 6, 0 )
500    ENDIF
501#endif
502
503    array_position  = 65536          !> Start offset for writing 2-D and 3.D arrays at 64 k
504    header_position = 0
505
506    header_int_index   = 1
507    header_char_index  = 1
508    header_real_index  = 1
509    header_array_index = 1
510
511    int_names    = ' '
512    int_values   = 0
513    text_lines   = ' '
514    real_names   = ' '
515    real_values  = 0.0
516    array_names  = ' '
517    array_offset = 0
518
519    int_names(1)     = 'nx'
520    int_values(1)    = nx
521    int_names(2)     = 'ny'
522    int_values(2)    = ny
523    int_names(3)     = 'nz'
524    int_values(3)    = nz
525    header_int_index = header_int_index+3
526
527    DO  i = 1, max_nr_arrays
528       array_offset(i) = 0
529       array_names(i)  = ' '
530    ENDDO
531
532    gh_size = STORAGE_SIZE( tgh ) / 8
533
534    IF ( rd_flag )  THEN
535       IF ( sm_io%iam_io_pe )  THEN
536!
537!--       File is open for reading
538#if defined( __parallel )
539!--       Set the default view
540          CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
541!
542!--       Read the file header size
543          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
544          CALL MPI_FILE_READ( fh, tgh, gh_size, MPI_BYTE, status, ierr )
545#else
546          CALL posix_lseek( fh, header_position )
547          buf_ptr = C_LOC( tgh )
548          CALL posix_read( fh, buf_ptr, gh_size )
549#endif
550       ENDIF
551#if defined( __parallel )
552       IF ( sm_io%is_sm_active() )  THEN
553          CALL MPI_BCAST( tgh, gh_size, MPI_BYTE, 0, sm_io%comm_shared, ierr )
554       ENDIF
555#endif
556       header_position = header_position + gh_size
557
558       include_total_domain_boundaries = ( tgh%i_outer_bound == 1 )
559
560!
561!--    File types depend on boundaries of the total domain being included in data. This has been
562!--    checked with the previous statement.
563       IF ( set_filetype )  THEN
564          CALL rd_mpi_io_create_filetypes
565          filetypes_created = .TRUE.
566       ENDIF
567
568       IF ( sm_io%iam_io_pe )  THEN
569#if defined( __parallel )
570!
571!--       Read INTEGER values
572          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
573          CALL MPI_FILE_READ( fh, int_names, SIZE( int_names ) * 32, MPI_CHAR, status, ierr )
574          header_position = header_position + SIZE( int_names ) * 32
575
576          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
577          CALL MPI_FILE_READ (fh, int_values, SIZE( int_values ), MPI_INT, status, ierr )
578          header_position = header_position + SIZE( int_values ) * iwp
579!
580!--       Character entries
581          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
582          CALL MPI_FILE_READ( fh, text_lines, SIZE( text_lines ) * 128, MPI_CHAR, status, ierr )
583          header_position = header_position + SIZE ( text_lines ) * 128
584!
585!--       REAL values
586          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
587          CALL MPI_FILE_READ( fh, real_names, SIZE( real_names ) * 32, MPI_CHAR, status, ierr )
588          header_position = header_position + SIZE( real_names ) * 32
589
590          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
591          CALL MPI_FILE_READ( fh, real_values, SIZE( real_values ), MPI_REAL, status, ierr )
592          header_position = header_position + SIZE( real_values ) * wp
593!
594!--       2d- and 3d-array headers
595          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
596          CALL MPI_FILE_READ( fh, array_names, SIZE( array_names ) * 32, MPI_CHAR, status, ierr )
597          header_position = header_position + SIZE( array_names ) * 32
598
599          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
600          CALL MPI_FILE_READ( fh, array_offset, SIZE( array_offset ) * MPI_OFFSET_KIND, MPI_BYTE,  &
601                              status,ierr )   ! there is no I*8 datatype in Fortran
602          header_position = header_position + SIZE( array_offset ) * rd_offset_kind
603#else
604          CALL posix_lseek( fh, header_position )
605          CALL posix_read( fh, int_names )
606          header_position = header_position + SIZE( int_names ) * 32
607
608          CALL posix_lseek( fh, header_position )
609          CALL posix_read( fh, int_values, SIZE( int_values ) )
610          header_position = header_position + SIZE( int_values ) * iwp
611!
612!--       Character entries
613          CALL posix_lseek( fh, header_position )
614          CALL posix_read( fh, text_lines )
615          header_position = header_position + SIZE( text_lines ) * 128
616!
617!--       REAL values
618          CALL posix_lseek( fh, header_position )
619          CALL posix_read( fh, real_names )
620          header_position = header_position + SIZE( real_names ) * 32
621
622          CALL posix_lseek( fh, header_position )
623          CALL posix_read( fh, real_values, SIZE( real_values ) )
624          header_position = header_position + SIZE( real_values ) * wp
625!
626!--       2d- and 3d-array headers
627          CALL posix_lseek( fh, header_position )
628          CALL posix_read( fh, array_names )
629          header_position = header_position + SIZE( array_names ) * 32
630
631          CALL posix_lseek( fh, header_position )
632          CALL posix_read( fh, array_offset, SIZE( array_offset ) ) ! there is no I*8 datatype in Fortran
633          header_position = header_position + SIZE( array_offset ) * rd_offset_kind
634#endif
635          IF ( debug_output )  CALL rd_mpi_io_print_header
636
637       ENDIF
638
639#if defined( __parallel )
640!
641!--    Broadcast header to all remaining cores that are not involved in I/O
642       IF ( sm_io%is_sm_active() )  THEN
643!
644!--        Not sure, that it is possible to broadcast CHARACTER array in one MPI_Bcast call
645           DO  i = 1, SIZE( int_names )
646              CALL MPI_BCAST( int_names(i), 32, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
647           ENDDO
648           CALL MPI_BCAST( int_values, SIZE( int_values ), MPI_INTEGER, 0, sm_io%comm_shared, ierr )
649
650           DO  i = 1, SIZE( text_lines )
651              CALL MPI_BCAST( text_lines(i), 128, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
652           ENDDO
653
654           DO  i = 1, SIZE( real_names )
655              CALL MPI_BCAST( real_names(i), 32, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
656           ENDDO
657           CALL MPI_BCAST( real_values, SIZE( real_values ), MPI_REAL, 0, sm_io%comm_shared, ierr )
658
659           DO  i = 1, SIZE( array_names )
660              CALL MPI_BCAST( array_names(i), 32, MPI_CHARACTER, 0, sm_io%comm_shared, ierr )
661           ENDDO
662           CALL MPI_BCAST( array_offset, SIZE( array_offset )*8, MPI_BYTE, 0, sm_io%comm_shared,   &
663                           ierr )  ! there is no I*8 datatype in Fortran (array_offset is I*8!)
664
665           CALL MPI_BCAST( header_position, rd_offset_kind, MPI_BYTE, 0, sm_io%comm_shared, ierr )
666
667       ENDIF
668#endif
669
670    ENDIF
671
672 END SUBROUTINE rd_mpi_io_open
673
674
675!--------------------------------------------------------------------------------------------------!
676! Description:
677! ------------
678!> Check, if array exists in restart file
679!--------------------------------------------------------------------------------------------------!
680 SUBROUTINE rd_mpi_io_check_array( name, found )
681
682    IMPLICIT NONE
683
684    CHARACTER(LEN=*), INTENT(IN) ::  name  !<
685
686    INTEGER(iwp) ::  i  !<
687
688    LOGICAl ::  found  !<
689
690
691    DO  i = 1, tgh%nr_arrays
692       IF ( TRIM( array_names(i) ) == TRIM( name ) )  THEN
693          array_position = array_offset(i)
694          found = .TRUE.
695          RETURN
696       ENDIF
697    ENDDO
698
699    found = .FALSE.
700
701 END SUBROUTINE rd_mpi_io_check_array
702
703
704
705!--------------------------------------------------------------------------------------------------!
706! Description:
707! ------------
708!> Read INTEGER with MPI-IO
709!--------------------------------------------------------------------------------------------------!
710 SUBROUTINE rrd_mpi_io_int( name, value )
711
712    IMPLICIT NONE
713
714    CHARACTER(LEN=*), INTENT(IN) :: name  !<
715
716    INTEGER(iwp)                   ::  i      !<
717    INTEGER(KIND=iwp), INTENT(OUT) ::  value  !<
718
719    LOGICAL ::  found  !<
720
721
722    found = .FALSE.
723    value = 0
724
725    DO  i = 1, tgh%nr_int
726       IF ( TRIM(int_names(i)) == TRIM( name ) )  THEN
727          value = int_values(i)
728          found = .TRUE.
729          EXIT
730       ENDIF
731    ENDDO
732
733    IF ( .NOT. found )  THEN
734       message_string = 'INTEGER variable "' // TRIM( name ) // '" not found in restart file'
735       CALL message( 'rrd_mpi_io_int', 'PA0722', 3, 2, 0, 6, 0 )
736    ENDIF
737
738 END SUBROUTINE rrd_mpi_io_int
739
740
741
742!--------------------------------------------------------------------------------------------------!
743! Description:
744! ------------
745!> Read REAL with MPI-IO
746!--------------------------------------------------------------------------------------------------!
747 SUBROUTINE rrd_mpi_io_real( name, value )
748
749    IMPLICIT NONE
750
751    CHARACTER(LEN=*), INTENT(IN) ::  name   !<
752
753    INTEGER(iwp)                 ::  i      !<
754
755    LOGICAL                      ::  found  !<
756
757    REAL(KIND=wp), INTENT(OUT)   ::  value  !<
758
759
760    found = .FALSE.
761    value = 0.0
762
763    DO  i = 1, tgh%nr_real
764       IF ( TRIM(real_names(i)) == TRIM( name ) )  THEN
765          value = real_values(i)
766          found = .TRUE.
767          EXIT
768       ENDIF
769    ENDDO
770
771    IF ( .NOT. found )  THEN
772       message_string = 'REAL variable "' // TRIM( name ) // '" not found in restart file'
773       CALL message( 'rrd_mpi_io_int', 'PA0722', 3, 2, 0, 6, 0 )
774    ENDIF
775
776 END SUBROUTINE rrd_mpi_io_real
777
778
779
780!--------------------------------------------------------------------------------------------------!
781! Description:
782! ------------
783!> Read 2d-real array with MPI-IO
784!--------------------------------------------------------------------------------------------------!
785 SUBROUTINE rrd_mpi_io_real_2d( name, data )
786
787    IMPLICIT NONE
788
789    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
790
791#if defined( __parallel )
792    INTEGER, DIMENSION(rd_status_size) ::  status  !<
793#endif
794    INTEGER(iwp)                       ::  i       !<
795
796    LOGICAL                            ::  found   !<
797
798    REAL(wp), INTENT(INOUT), DIMENSION(nysg:nyng,nxlg:nxrg) ::  data  !<
799
800
801    found = .FALSE.
802
803    DO  i = 1, tgh%nr_arrays
804       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
805          array_position = array_offset(i)
806          found = .TRUE.
807          EXIT
808       ENDIF
809    ENDDO
810
811    IF ( found )  THEN
812#if defined( __parallel )
813       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
814       IF ( sm_io%iam_io_pe )  THEN
815          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL,    &
816                                  ierr )
817          CALL MPI_FILE_READ_ALL( fh, array_2d, SIZE( array_2d ), MPI_REAL, status, ierr )
818       ENDIF
819       CALL sm_io%sm_node_barrier()
820#else
821       CALL posix_lseek( fh, array_position )
822       CALL posix_read( fh, array_2d, SIZE( array_2d ) )
823#endif
824
825       IF ( include_total_domain_boundaries)  THEN
826          DO  i = lb%nxl, lb%nxr
827             data(lb%nys-nbgp:lb%nyn-nbgp,i-nbgp) = array_2d(i,lb%nys:lb%nyn)
828          ENDDO
829          IF ( debug_level >= 2)  WRITE(9,*) 'r2f_ob ', TRIM(name),' ', SUM( data(nys:nyn,nxl:nxr) )
830       ELSE
831          DO  i = nxl, nxr
832             data(nys:nyn,i) = array_2d(i,nys:nyn)
833          ENDDO
834          IF ( debug_level >= 2) WRITE(9,*) 'r2f ', TRIM( name ),' ', SUM( data(nys:nyn,nxl:nxr) )
835       ENDIF
836
837       CALL exchange_horiz_2d( data )
838
839    ELSE
840       message_string = '2d-REAL array "' // TRIM( name ) // '" not found in restart file'
841       CALL message( 'rrd_mpi_io_int', 'PA0722', 3, 2, 0, 6, 0 )
842    ENDIF
843
844 END SUBROUTINE rrd_mpi_io_real_2d
845
846
847
848!--------------------------------------------------------------------------------------------------!
849! Description:
850! ------------
851!> Read 2d-INTEGER array with MPI-IO
852!--------------------------------------------------------------------------------------------------!
853 SUBROUTINE rrd_mpi_io_int_2d( name, data )
854
855    IMPLICIT NONE
856
857    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
858
859    INTEGER(iwp)                       ::  i       !<
860    INTEGER(iwp)                       ::  j       !<
861
862#if defined( __parallel )
863    INTEGER, DIMENSION(rd_status_size) ::  status  !<
864#endif
865
866    INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:,:) ::  data  !<
867
868    LOGICAL ::  found  !<
869
870
871    found = .FALSE.
872
873    DO  i = 1, tgh%nr_arrays
874       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
875          array_position = array_offset(i)
876          found = .TRUE.
877          EXIT
878       ENDIF
879    ENDDO
880
881    IF ( found )  THEN
882
883       IF ( ( nxr - nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) )  THEN
884!
885!--       Output array with Halos.
886!--       ATTENTION: INTEGER arrays with ghost boundaries are not implemented yet. This kind of
887!--                  array would be dimensioned in the caller subroutine like this:
888!--                  INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg)::  data
889          message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be read from restart ' // &
890                           'file is defined with illegal dimensions in the PALM code'
891          CALL message( 'rrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
892
893       ELSEIF ( (nxr - nxl + 1) == SIZE( data, 2 ) )  THEN
894!
895!--       INTEGER input array without Halos.
896!--       This kind of array is dimensioned in the caller subroutine
897!--       INTEGER, DIMENSION(nys:nyn,nxl:nxr) ::  data
898
899#if defined( __parallel )
900          CALL sm_io%sm_node_barrier() ! Has no effect if I/O on limited number of cores is inactive
901          IF ( sm_io%iam_io_pe )  THEN
902             CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native',         &
903                                     MPI_INFO_NULL, ierr )
904             CALL MPI_FILE_READ_ALL( fh, array_2di, SIZE( array_2di ), MPI_INTEGER, status, ierr )
905          ENDIF
906          CALL sm_io%sm_node_barrier()
907#else
908          CALL posix_lseek( fh, array_position )
909          CALL posix_read( fh, array_2di, SIZE( array_2di ) )
910#endif
911
912          DO  j = nys, nyn
913             DO  i = nxl, nxr
914                data(j-nys+1,i-nxl+1) = array_2di(i,j)
915             ENDDO
916          ENDDO
917
918       ELSE
919
920          message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be read from restart ' // &
921                           'file is defined with illegal dimensions in the PALM code'
922          CALL message( 'rrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
923
924       ENDIF
925
926    ELSE
927
928       message_string = '2d-INTEGER array "' // TRIM( name ) // '" not found in restart file'
929       CALL message( 'rrd_mpi_io_int_2d', 'PA0722', 3, 2, 0, 6, 0 )
930
931    ENDIF
932
933 END SUBROUTINE rrd_mpi_io_int_2d
934
935
936
937!--------------------------------------------------------------------------------------------------!
938! Description:
939! ------------
940!> Read 2d-REAL array with MPI-IO
941!--------------------------------------------------------------------------------------------------!
942 SUBROUTINE rrd_mpi_io_real_3d( name, data )
943
944    IMPLICIT NONE
945
946    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
947
948    INTEGER(iwp)                       ::  i       !<
949
950#if defined( __parallel )
951    INTEGER, DIMENSION(rd_status_size) ::  status  !<
952#endif
953
954    LOGICAL                            ::  found   !<
955
956    REAL(wp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data  !<
957
958
959    found = .FALSE.
960
961    DO  i = 1, tgh%nr_arrays
962       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
963          array_position = array_offset(i)
964          found = .TRUE.
965          EXIT
966       ENDIF
967    ENDDO
968
969    IF ( found )  THEN
970#if defined( __parallel )
971       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
972       IF( sm_io%iam_io_pe )  THEN
973          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL,    &
974                                  ierr )
975          CALL MPI_FILE_READ_ALL( fh, array_3d, SIZE( array_3d ), MPI_REAL, status, ierr )
976       ENDIF
977       CALL sm_io%sm_node_barrier()
978#else
979       CALL posix_lseek( fh, array_position )
980       CALL posix_read(fh, array_3d, SIZE( array_3d ) )
981#endif
982       IF ( include_total_domain_boundaries)  THEN
983          DO  i = lb%nxl, lb%nxr
984             data(:,lb%nys-nbgp:lb%nyn-nbgp,i-nbgp) = array_3d(:,i,lb%nys:lb%nyn)
985          ENDDO
986       ELSE
987          DO  i = nxl, nxr
988             data(:,nys:nyn,i) = array_3d(:,i,nys:nyn)
989          ENDDO
990       ENDIF
991
992       CALL exchange_horiz( data, nbgp )
993
994    ELSE
995
996       message_string = '3d-REAL array "' // TRIM( name ) // '" not found in restart file'
997       CALL message( 'rrd_mpi_io_real_3d', 'PA0722', 3, 2, 0, 6, 0 )
998
999    ENDIF
1000
1001 END SUBROUTINE rrd_mpi_io_real_3d
1002
1003
1004
1005!--------------------------------------------------------------------------------------------------!
1006! Description:
1007! ------------
1008!> Read 3d-REAL soil array with MPI-IO
1009!> nzb_soil, nzt_soil are located in the module land_surface_model_mod. Since Fortran does not allow
1010!> cross referencing of module variables, it is required to pass these variables as arguments.
1011!--------------------------------------------------------------------------------------------------!
1012 SUBROUTINE rrd_mpi_io_real_3d_soil( name, data, nzb_soil, nzt_soil )
1013
1014    IMPLICIT NONE
1015
1016    CHARACTER(LEN=*), INTENT(IN)       ::  name      !<
1017
1018    INTEGER(iwp)                       ::  i         !<
1019    INTEGER, INTENT(IN)                ::  nzb_soil  !<
1020    INTEGER, INTENT(IN)                ::  nzt_soil  !<
1021
1022#if defined( __parallel )
1023    INTEGER, DIMENSION(rd_status_size) ::  status    !<
1024#endif
1025
1026    LOGICAL                            ::  found     !<
1027
1028    REAL(wp), INTENT(INOUT), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ::  data  !<
1029
1030
1031    found = .FALSE.
1032
1033    DO  i = 1, tgh%nr_arrays
1034       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
1035          array_position = array_offset(i)
1036          found = .TRUE.
1037          EXIT
1038       ENDIF
1039    ENDDO
1040
1041    IF ( found )  THEN
1042#if defined( __parallel )
1043       CALL rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
1044       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
1045       IF ( sm_io%iam_io_pe )  THEN
1046          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native',               &
1047                                  MPI_INFO_NULL, ierr )
1048          CALL MPI_FILE_READ_ALL( fh, array_3d_soil, SIZE( array_3d_soil ), MPI_REAL, status, ierr )
1049          CALL MPI_TYPE_FREE( ft_3dsoil, ierr )
1050       ENDIF
1051       CALL sm_io%sm_node_barrier()
1052#else
1053       CALL posix_lseek( fh, array_position )
1054       CALL posix_read( fh, array_3d_soil, SIZE( array_3d_soil ) )
1055#endif
1056       IF ( include_total_domain_boundaries )  THEN
1057          DO  i = lb%nxl, lb%nxr
1058             data(:,lb%nys-nbgp:lb%nyn-nbgp,i-nbgp) = array_3d(:,i,lb%nys:lb%nyn)
1059          ENDDO
1060       ELSE
1061          DO  i = nxl, nxr
1062             data(:,nys:nyn,i) = array_3d(:,i,nys:nyn)
1063          ENDDO
1064       ENDIF
1065
1066    ELSE
1067
1068       message_string = '3d-REAL soil array "' // TRIM( name ) // '" not found in restart file'
1069       CALL message( 'rrd_mpi_io_real_3d_soil', 'PA0722', 3, 2, 0, 6, 0 )
1070
1071    ENDIF
1072
1073 END SUBROUTINE rrd_mpi_io_real_3d_soil
1074
1075
1076
1077!--------------------------------------------------------------------------------------------------!
1078! Description:
1079! ------------
1080!> Read CHARACTER with MPI-IO
1081!--------------------------------------------------------------------------------------------------!
1082 SUBROUTINE rrd_mpi_io_char( name, text )
1083
1084    IMPLICIT NONE
1085
1086    CHARACTER(LEN=*), INTENT(IN)  ::  name   !<
1087    CHARACTER(LEN=*), INTENT(OUT) ::  text   !<
1088    CHARACTER(LEN=128)            ::  line   !<
1089
1090    INTEGER(iwp)                  ::  i      !<
1091
1092    LOGICAL                       ::  found  !<
1093
1094
1095    found = .FALSE.
1096    text = ' '
1097
1098    DO  i = 1, tgh%nr_char
1099       line = text_lines(i)
1100       IF ( TRIM( line(1:32) ) == TRIM( name ) )  THEN
1101          text = line(33:)
1102          found = .TRUE.
1103          EXIT
1104       ENDIF
1105    ENDDO
1106
1107    IF ( .NOT. found )  THEN
1108       message_string = 'CHARACTER variable "' // TRIM( name ) // '" not found in restart file'
1109       CALL message( 'rrd_mpi_io_char', 'PA0722', 3, 2, 0, 6, 0 )
1110    ENDIF
1111
1112 END SUBROUTINE rrd_mpi_io_char
1113
1114
1115
1116!--------------------------------------------------------------------------------------------------!
1117! Description:
1118! ------------
1119!> Read LOGICAL with MPI-IO
1120!--------------------------------------------------------------------------------------------------!
1121 SUBROUTINE rrd_mpi_io_logical( name, value )
1122
1123    IMPLICIT NONE
1124
1125    CHARACTER(LEN=*), INTENT(IN) ::  name                !<
1126
1127    INTEGER(iwp)                 ::  logical_as_integer  !<
1128
1129    LOGICAL, INTENT(OUT)         ::  value               !<
1130
1131
1132    CALL rrd_mpi_io_int( name, logical_as_integer )
1133    value = ( logical_as_integer == 1 )
1134
1135 END SUBROUTINE rrd_mpi_io_logical
1136
1137
1138
1139!--------------------------------------------------------------------------------------------------!
1140! Description:
1141! ------------
1142!> Write INTEGER with MPI-IO
1143!--------------------------------------------------------------------------------------------------!
1144 SUBROUTINE wrd_mpi_io_int( name, value )
1145
1146    IMPLICIT NONE
1147
1148    CHARACTER(LEN=*), INTENT(IN)  ::  name   !<
1149
1150    INTEGER(KIND=iwp), INTENT(IN) ::  value  !<
1151
1152
1153    IF ( header_int_index == max_nr_int )  THEN
1154       STOP '+++ maximum number of INTEGER entries in restart file header exceeded'
1155    ENDIF
1156
1157    int_names(header_int_index)  = name
1158    int_values(header_int_index) = value
1159    header_int_index = header_int_index + 1
1160
1161 END SUBROUTINE wrd_mpi_io_int
1162
1163
1164!--------------------------------------------------------------------------------------------------!
1165! Description:
1166! ------------
1167!> To do: Description missing!
1168!--------------------------------------------------------------------------------------------------!
1169 SUBROUTINE wrd_mpi_io_real( name, value )
1170
1171    IMPLICIT NONE
1172
1173    CHARACTER(LEN=*), INTENT(IN) ::  name   !<
1174
1175    REAL(wp), INTENT(IN)         ::  value  !<
1176
1177
1178    IF ( header_real_index == max_nr_real )  THEN
1179       STOP '+++ maximum number of REAL entries in restart file header exceeded'
1180    ENDIF
1181
1182    real_names(header_real_index)  = name
1183    real_values(header_real_index) = value
1184    header_real_index = header_real_index + 1
1185
1186 END SUBROUTINE wrd_mpi_io_real
1187
1188
1189
1190!--------------------------------------------------------------------------------------------------!
1191! Description:
1192! ------------
1193!> Write 2d-REAL array with MPI-IO
1194!--------------------------------------------------------------------------------------------------!
1195 SUBROUTINE wrd_mpi_io_real_2d( name, data )
1196
1197    IMPLICIT NONE
1198
1199    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
1200
1201    INTEGER(iwp)                       ::  i       !<
1202
1203#if defined( __parallel )
1204    INTEGER, DIMENSION(rd_status_size) ::  status  !<
1205#endif
1206
1207    REAL(wp), INTENT(IN), DIMENSION(nysg:nyng,nxlg:nxrg) ::  data  !<
1208
1209
1210    IF ( header_array_index == max_nr_arrays )  THEN
1211       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
1212    ENDIF
1213
1214    array_names(header_array_index)  = name
1215    array_offset(header_array_index) = array_position
1216    header_array_index = header_array_index + 1
1217
1218    IF ( include_total_domain_boundaries )  THEN
1219!
1220!--    Prepare output with outer boundaries
1221       DO  i = lb%nxl, lb%nxr
1222          array_2d(i,lb%nys:lb%nyn) = data(lb%nys-nbgp:lb%nyn-nbgp,i-nbgp)
1223       ENDDO
1224
1225    ELSE
1226!
1227!--    Prepare output without outer boundaries
1228       DO  i = nxl,nxr
1229          array_2d(i,lb%nys:lb%nyn) = data(nys:nyn,i)
1230       ENDDO
1231
1232    ENDIF
1233
1234#if defined( __parallel )
1235    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
1236    IF ( sm_io%iam_io_pe )  THEN
1237       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_2d, 'native', MPI_INFO_NULL, ierr )
1238       CALL MPI_FILE_WRITE_ALL( fh, array_2d, SIZE( array_2d), MPI_REAL, status, ierr )
1239    ENDIF
1240    CALL sm_io%sm_node_barrier()
1241#else
1242    CALL posix_lseek( fh, array_position )
1243    CALL posix_write( fh, array_2d, SIZE( array_2d ) )
1244#endif
1245!
1246!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
1247!-- Maybe a compiler problem.
1248    array_position = array_position + ( INT( lb%ny, KIND=rd_offset_kind ) + 1 ) *                  &
1249                                      ( INT( lb%nx, KIND=rd_offset_kind ) + 1 ) * wp
1250
1251 END SUBROUTINE wrd_mpi_io_real_2d
1252
1253
1254
1255!--------------------------------------------------------------------------------------------------!
1256! Description:
1257! ------------
1258!> Write 2d-INTEGER array with MPI-IO
1259!--------------------------------------------------------------------------------------------------!
1260 SUBROUTINE wrd_mpi_io_int_2d( name, data )
1261
1262    IMPLICIT NONE
1263
1264    CHARACTER(LEN=*), INTENT(IN)                  ::  name    !<
1265
1266    INTEGER(iwp)                                  ::  i       !<
1267    INTEGER(iwp)                                  ::  j       !<
1268
1269#if defined( __parallel )
1270    INTEGER, DIMENSION(rd_status_size)            ::  status  !<
1271#endif
1272    INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:,:) ::  data    !<
1273
1274
1275    IF ( header_array_index == max_nr_arrays )  THEN
1276       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
1277    ENDIF
1278
1279    array_names(header_array_index)  = name
1280    array_offset(header_array_index) = array_position
1281    header_array_index = header_array_index + 1
1282
1283    IF ( ( nxr - nxl + 1 + 2 * nbgp ) == SIZE( data, 2 ) )  THEN
1284!
1285!--    Integer arrays with ghost layers are not implemented yet. These kind of arrays would be
1286!--    dimensioned in the caller subroutine as
1287!--    INTEGER, DIMENSION(nysg:nyng,nxlg:nxrg) ::  data
1288       message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be written to restart ' //   &
1289                        'file is defined with illegal dimensions in the PALM code'
1290       CALL message( 'wrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
1291
1292    ELSEIF ( ( nxr-nxl+1 ) == SIZE( data, 2 ) )  THEN
1293!
1294!--    INTEGER input array without ghost layers.
1295!--    This kind of array is dimensioned in the caller subroutine as
1296!--    INTEGER, DIMENSION(nys:nyn,nxl:nxr) ::  data
1297       DO  j = nys, nyn
1298          DO  i = nxl, nxr
1299             array_2di(i,j) = data(j-nys+1,i-nxl+1)
1300          ENDDO
1301       ENDDO
1302#if defined( __parallel )
1303       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
1304       IF ( sm_io%iam_io_pe )  THEN
1305          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_2di_nb, 'native',            &
1306                                  MPI_INFO_NULL, ierr )
1307          CALL MPI_FILE_WRITE_ALL( fh, array_2di, SIZE( array_2di ), MPI_INTEGER, status, ierr )
1308       ENDIF
1309       CALL sm_io%sm_node_barrier()
1310#else
1311       CALL posix_lseek( fh, array_position )
1312       CALL posix_write( fh, array_2di, SIZE( array_2di ) )
1313#endif
1314!
1315!--    Type conversion required, otherwise rigth hand side brackets are calculated assuming 4 byte
1316!--    INT. Maybe a compiler problem.
1317       array_position = array_position + INT( (ny+1), KIND=rd_offset_kind ) *                      &
1318                                         INT( (nx+1), KIND=rd_offset_kind ) * 4
1319
1320    ELSE
1321
1322       message_string = '2d-INTEGER array "' // TRIM( name ) // '" to be written to restart ' //   &
1323                        'file is defined with illegal dimensions in the PALM code'
1324       CALL message( 'wrd_mpi_io_int_2d', 'PA0723', 3, 2, 0, 6, 0 )
1325
1326    ENDIF
1327
1328 END SUBROUTINE wrd_mpi_io_int_2d
1329
1330
1331
1332!--------------------------------------------------------------------------------------------------!
1333! Description:
1334! ------------
1335!> Write 3d-REAL array with MPI-IO
1336!--------------------------------------------------------------------------------------------------!
1337 SUBROUTINE wrd_mpi_io_real_3d( name, data )
1338
1339    IMPLICIT NONE
1340
1341    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
1342
1343    INTEGER(iwp)                       ::  i       !<
1344#if defined( __parallel )
1345    INTEGER, DIMENSION(rd_status_size) ::  status  !<
1346#endif
1347    REAL(wp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data !<
1348
1349
1350    IF ( header_array_index == max_nr_arrays )  THEN
1351       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
1352    ENDIF
1353
1354    array_names(header_array_index)  = name
1355    array_offset(header_array_index) = array_position
1356    header_array_index = header_array_index + 1
1357
1358    IF ( include_total_domain_boundaries )  THEN
1359!
1360!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
1361!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
1362!--    index order of the array in the same way, i.e. the first dimension should be along x and the
1363!--    second along y. For this reason, the original PALM data need to be swaped.
1364       DO  i = lb%nxl, lb%nxr
1365          array_3d(:,i,lb%nys:lb%nyn) = data(:,lb%nys-nbgp:lb%nyn-nbgp,i-nbgp)
1366       ENDDO
1367
1368    ELSE
1369!
1370!--    Prepare output of 3d-REAL-array without ghost layers
1371       DO  i = nxl, nxr
1372           array_3d(:,i,lb%nys:lb%nyn) = data(:,nys:nyn,i)
1373       ENDDO
1374
1375    ENDIF
1376#if defined( __parallel )
1377    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
1378    IF ( sm_io%iam_io_pe )  THEN
1379       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3d, 'native', MPI_INFO_NULL, ierr )
1380       CALL MPI_FILE_WRITE_ALL( fh, array_3d, SIZE( array_3d ), MPI_REAL, status, ierr )
1381    ENDIF
1382    CALL sm_io%sm_node_barrier()
1383#else
1384    CALL posix_lseek( fh, array_position )
1385    CALL posix_write( fh, array_3d, SIZE( array_3d ) )
1386#endif
1387!
1388!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
1389!-- Maybe a compiler problem.
1390    array_position = array_position + INT(    (nz+2), KIND = rd_offset_kind ) *                    &
1391                                      INT( (lb%ny+1), KIND = rd_offset_kind ) *                    &
1392                                      INT( (lb%nx+1), KIND = rd_offset_kind ) * wp
1393
1394 END SUBROUTINE wrd_mpi_io_real_3d
1395
1396
1397
1398!--------------------------------------------------------------------------------------------------!
1399! Description:
1400! ------------
1401!> Write 3d-REAL soil array with MPI-IO.
1402!> nzb_soil, nzt_soil are located in the module land_surface_model_mod. Since Fortran does not allow
1403!> cross referencing of module variables, it is required to pass these variables as arguments.
1404!--------------------------------------------------------------------------------------------------!
1405 SUBROUTINE wrd_mpi_io_real_3d_soil( name, data, nzb_soil, nzt_soil )
1406
1407    IMPLICIT NONE
1408
1409    CHARACTER(LEN=*), INTENT(IN)       ::  name      !<
1410
1411    INTEGER(iwp)                       ::  i         !<
1412    INTEGER, INTENT(IN)                ::  nzb_soil  !<
1413    INTEGER, INTENT(IN)                ::  nzt_soil  !<
1414
1415#if defined( __parallel )
1416    INTEGER, DIMENSION(rd_status_size) ::  status  !<
1417#endif
1418
1419    REAL(wp), INTENT(IN), DIMENSION(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) ::  data  !<
1420
1421
1422    IF ( header_array_index == max_nr_arrays )  THEN
1423       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
1424    ENDIF
1425
1426    array_names(header_array_index)  = name
1427    array_offset(header_array_index) = array_position
1428    header_array_index = header_array_index + 1
1429
1430#if defined( __parallel )
1431    CALL rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
1432#endif
1433
1434    IF ( include_total_domain_boundaries)  THEN
1435!
1436!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
1437!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
1438!--    index order of the array in the same way, i.e. the first dimension should be along x and the
1439!--    second along y. For this reason, the original PALM data need to be swaped.
1440       DO  i = lb%nxl, lb%nxr
1441          array_3d_soil(:,i,lb%nys:lb%nyn) = data(:,lb%nys-nbgp:lb%nyn-nbgp,i-nbgp)
1442       ENDDO
1443
1444    ELSE
1445!
1446!--    Prepare output of 3d-REAL-array without ghost layers
1447       DO  i = nxl, nxr
1448          array_3d_soil(:,i,lb%nys:lb%nyn) = data(:,nys:nyn,i)
1449       ENDDO
1450
1451    ENDIF
1452#if defined( __parallel )
1453    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
1454    IF ( sm_io%iam_io_pe )  THEN
1455       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_3dsoil, 'native', MPI_INFO_NULL,   &
1456                               ierr )
1457       CALL MPI_FILE_WRITE_ALL( fh, array_3d_soil, SIZE( array_3d_soil ), MPI_REAL, status, ierr )
1458    ENDIF
1459    CALL sm_io%sm_node_barrier()
1460#else
1461    CALL posix_lseek( fh, array_position )
1462    CALL posix_write( fh, array_3d_soil, SIZE( array_3d_soil ) )
1463#endif
1464!
1465!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
1466!-- Maybe a compiler problem.
1467    array_position = array_position + INT( (nzt_soil-nzb_soil+1), KIND = rd_offset_kind ) *        &
1468                                      INT( (lb%ny+1),             KIND = rd_offset_kind ) *        &
1469                                      INT( (lb%nx+1),             KIND = rd_offset_kind ) * wp
1470
1471 END SUBROUTINE wrd_mpi_io_real_3d_soil
1472
1473
1474
1475!--------------------------------------------------------------------------------------------------!
1476! Description:
1477! ------------
1478!> Write CHARATCTER with MPI-IO
1479!--------------------------------------------------------------------------------------------------!
1480 SUBROUTINE wrd_mpi_io_char( name, text )
1481
1482    IMPLICIT NONE
1483
1484    CHARACTER(LEN=128)           ::  lo_line  !<
1485    CHARACTER(LEN=*), INTENT(IN) ::  name     !<
1486    CHARACTER(LEN=*), INTENT(IN) ::  text     !<
1487
1488
1489    IF ( header_char_index == max_nr_char )  THEN
1490       STOP '+++ maximum number of CHARACTER entries in restart file header exceeded'
1491    ENDIF
1492
1493    lo_line      = name
1494    lo_line(33:) = text
1495    text_lines(header_char_index) = lo_line
1496    header_char_index = header_char_index + 1
1497
1498 END SUBROUTINE wrd_mpi_io_char
1499
1500
1501
1502!--------------------------------------------------------------------------------------------------!
1503! Description:
1504! ------------
1505!> Write LOGICAL with MPI-IO
1506!--------------------------------------------------------------------------------------------------!
1507 SUBROUTINE wrd_mpi_io_logical( name, value )
1508
1509    IMPLICIT NONE
1510
1511    CHARACTER(LEN=*), INTENT(IN) ::  name                !<
1512
1513    INTEGER(iwp)                 ::  logical_as_integer  !<
1514
1515    LOGICAL, INTENT(IN)          ::  value               !<
1516
1517
1518    IF ( value )  THEN
1519       logical_as_integer = 1
1520    ELSE
1521       logical_as_integer = 0
1522    ENDIF
1523
1524    CALL wrd_mpi_io_int( name, logical_as_integer )
1525
1526 END SUBROUTINE wrd_mpi_io_logical
1527
1528
1529
1530!--------------------------------------------------------------------------------------------------!
1531! Description:
1532! ------------
1533!> Read 1d-REAL global array with MPI-IO.
1534!> Array contains identical data on all PEs.
1535!--------------------------------------------------------------------------------------------------!
1536 SUBROUTINE rrd_mpi_io_global_array_real_1d( name, data )
1537
1538    IMPLICIT NONE
1539
1540    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
1541
1542    INTEGER(iwp)                       ::  i       !<
1543    INTEGER(KIND=rd_offset_kind)       ::  offset  !<
1544
1545#if defined( __parallel )
1546    INTEGER, DIMENSION(rd_status_size) ::  status  !<
1547#endif
1548
1549    LOGICAL                            ::  found   !<
1550
1551    REAL(KIND=wp), INTENT(INOUT), DIMENSION(:) ::  data  !<
1552
1553
1554    offset = 0
1555    found  = .FALSE.
1556
1557    DO  i = 1, tgh%nr_arrays
1558       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
1559          array_position = array_offset(i)
1560          found = .TRUE.
1561          EXIT
1562       ENDIF
1563    ENDDO
1564
1565    IF ( found )  THEN
1566!
1567!--    Set default view
1568#if defined( __parallel )
1569       IF ( sm_io%iam_io_pe )  THEN
1570          CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
1571          CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
1572          CALL MPI_FILE_READ_ALL( fh, data, SIZE( data ), MPI_REAL, status, ierr )
1573       ENDIF
1574       IF ( sm_io%is_sm_active() )  THEN
1575          CALL MPI_BCAST( data, SIZE( data ), MPI_REAL, 0, sm_io%comm_shared, ierr )
1576       ENDIF
1577#else
1578       CALL posix_lseek( fh, array_position )
1579       CALL posix_read( fh, data, SIZE( data ) )
1580#endif
1581
1582    ELSE
1583
1584       message_string = '1d/2d/3d/4d-REAL global array "' // TRIM( name ) // '" not found in ' //  &
1585                        'restart file'
1586       CALL message( 'rrd_mpi_io_global_array_real_1d', 'PA0722', 3, 2, 0, 6, 0 )
1587
1588    ENDIF
1589
1590 END SUBROUTINE rrd_mpi_io_global_array_real_1d
1591
1592
1593
1594!--------------------------------------------------------------------------------------------------!
1595! Description:
1596! ------------
1597!> Read 2d-REAL global array with MPI-IO.
1598!> Array contains identical data on all PEs.
1599!--------------------------------------------------------------------------------------------------!
1600 SUBROUTINE rrd_mpi_io_global_array_real_2d( name, data )
1601
1602    IMPLICIT NONE
1603
1604    CHARACTER(LEN=*), INTENT(IN)                      ::  name      !<
1605
1606    INTEGER, DIMENSION(1)                             ::  bufshape  !<
1607
1608    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET ::  data      !<
1609    REAL(KIND=wp), POINTER, DIMENSION(:)              ::  buf       !<
1610
1611    TYPE(C_PTR)                                       ::  c_data    !<
1612
1613
1614    c_data = C_LOC( data )
1615    bufshape(1) = SIZE( data )
1616    CALL C_F_POINTER( c_data, buf, bufshape )
1617
1618    CALL rrd_mpi_io_global_array_real_1d( name, buf )
1619
1620 END SUBROUTINE rrd_mpi_io_global_array_real_2d
1621
1622
1623
1624!--------------------------------------------------------------------------------------------------!
1625! Description:
1626! ------------
1627!> Read 3d-REAL global array with MPI-IO.
1628!> Array contains identical data on all PEs.
1629!--------------------------------------------------------------------------------------------------!
1630 SUBROUTINE rrd_mpi_io_global_array_real_3d( name, data )
1631
1632    IMPLICIT NONE
1633
1634    CHARACTER(LEN=*), INTENT(IN)                        ::  name      !<
1635
1636    INTEGER, DIMENSION(1)                               ::  bufshape  !<
1637
1638    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET ::  data      !<
1639    REAL(KIND=wp), POINTER, DIMENSION(:)                ::  buf       !<
1640
1641    TYPE(C_PTR)                                         ::  c_data    !<
1642
1643
1644    c_data = C_LOC( data )
1645    bufshape(1) = SIZE( data )
1646    CALL C_F_POINTER( c_data, buf, bufshape )
1647
1648    CALL rrd_mpi_io_global_array_real_1d( name, buf )
1649
1650 END SUBROUTINE rrd_mpi_io_global_array_real_3d
1651
1652
1653
1654!--------------------------------------------------------------------------------------------------!
1655! Description:
1656! ------------
1657!> Read 4d-REAL global array with MPI-IO.
1658!> Array contains identical data on all PEs.
1659!--------------------------------------------------------------------------------------------------!
1660 SUBROUTINE rrd_mpi_io_global_array_real_4d( name, data )
1661
1662    IMPLICIT NONE
1663
1664    CHARACTER(LEN=*), INTENT(IN)                          ::  name      !<
1665
1666    INTEGER, DIMENSION(1)                                 ::  bufshape  !<
1667
1668    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET ::  data      !<
1669    REAL(KIND=wp), POINTER, DIMENSION(:)                  ::  buf       !<
1670
1671    TYPE(C_PTR)                                           ::  c_data    !<
1672
1673
1674    c_data = C_LOC( data )
1675    bufshape(1) = SIZE( data)
1676    CALL C_F_POINTER( c_data, buf, bufshape )
1677
1678    CALL rrd_mpi_io_global_array_real_1d( name, buf )
1679
1680 END SUBROUTINE rrd_mpi_io_global_array_real_4d
1681
1682
1683
1684!--------------------------------------------------------------------------------------------------!
1685! Description:
1686! ------------
1687!> Read 1d-INTEGER global array with MPI-IO.
1688!> Array contains identical data on all PEs.
1689!--------------------------------------------------------------------------------------------------!
1690 SUBROUTINE rrd_mpi_io_global_array_int_1d( name, data )
1691
1692    IMPLICIT NONE
1693
1694    CHARACTER(LEN=*), INTENT(IN)                   ::  name    !<
1695
1696    INTEGER(iwp)                                   ::  i       !<
1697    INTEGER(KIND=rd_offset_kind)                   ::  offset  !<
1698
1699#if defined( __parallel )
1700    INTEGER, DIMENSION(rd_status_size)             ::  status  !<
1701#endif
1702    INTEGER(KIND=iwp), INTENT(INOUT), DIMENSION(:) ::  data    !<
1703
1704    LOGICAL                                        ::  found   !<
1705
1706
1707    offset = 0
1708    found  = .FALSE.
1709
1710    DO  i = 1, tgh%nr_arrays
1711       IF ( TRIM( array_names(i) ) == TRIM( name ) )  THEN
1712          array_position = array_offset(i)
1713          found = .TRUE.
1714          EXIT
1715       ENDIF
1716    ENDDO
1717
1718    IF ( found )  THEN
1719!
1720!--    Set default view
1721#if defined( __parallel )
1722       IF ( sm_io%iam_io_pe )  THEN
1723          CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
1724          CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
1725          CALL MPI_FILE_READ_ALL( fh, data, SIZE( data), MPI_INTEGER, status, ierr )
1726       ENDIF
1727       IF ( sm_io%is_sm_active() )  THEN
1728          CALL MPI_BCAST( data, SIZE( data ), MPI_INTEGER, 0, sm_io%comm_shared, ierr )
1729       ENDIF
1730#else
1731       CALL posix_lseek( fh, array_position )
1732       CALL posix_read( fh, data, SIZE( data ) )
1733#endif
1734    ELSE
1735
1736       message_string = '1d-INTEGER global array "' // TRIM( name ) // '" not found in ' //        &
1737                        'restart file'
1738       CALL message( 'rrd_mpi_io_global_array_int_1d', 'PA0722', 3, 2, 0, 6, 0 )
1739
1740    ENDIF
1741
1742 END SUBROUTINE rrd_mpi_io_global_array_int_1d
1743
1744
1745
1746!--------------------------------------------------------------------------------------------------!
1747! Description:
1748! ------------
1749!> Write 1d-REAL global array with MPI-IO.
1750!> Array contains identical data on all PEs.
1751!--------------------------------------------------------------------------------------------------!
1752 SUBROUTINE wrd_mpi_io_global_array_real_1d( name, data )
1753
1754    IMPLICIT NONE
1755
1756    CHARACTER(LEN=*), INTENT(IN)            ::  name    !<
1757
1758    INTEGER(KIND=rd_offset_kind)            ::  offset  !<
1759
1760#if defined( __parallel )
1761    INTEGER, DIMENSION(rd_status_size)      ::  status  !<
1762#endif
1763
1764    REAL(KIND=wp), INTENT(IN), DIMENSION(:) ::  data    !<
1765
1766
1767    offset = 0
1768
1769    IF ( header_array_index == max_nr_arrays )  THEN
1770       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
1771    ENDIF
1772
1773    array_names(header_array_index)  = name
1774    array_offset(header_array_index) = array_position
1775    header_array_index = header_array_index + 1
1776
1777!
1778!-- Set default view
1779#if defined( __parallel )
1780    IF ( sm_io%iam_io_pe )  THEN
1781       CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
1782    ENDIF
1783!
1784!-- Only PE 0 writes replicated data
1785    IF ( myid == 0 )  THEN
1786       CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
1787       CALL MPI_FILE_WRITE( fh, data, SIZE( data ), MPI_REAL, status, ierr )
1788    ENDIF
1789#else
1790    CALL posix_lseek( fh, array_position )
1791    CALL posix_write( fh, data, SIZE( data ) )
1792#endif
1793    array_position = array_position + SIZE( data ) * wp
1794
1795 END SUBROUTINE wrd_mpi_io_global_array_real_1d
1796
1797
1798
1799!--------------------------------------------------------------------------------------------------!
1800! Description:
1801! ------------
1802!> Write 2d-REAL global array with MPI-IO.
1803!> Array contains identical data on all PEs.
1804!--------------------------------------------------------------------------------------------------!
1805 SUBROUTINE wrd_mpi_io_global_array_real_2d( name, data )
1806
1807    IMPLICIT NONE
1808
1809    CHARACTER(LEN=*), INTENT(IN)                      ::  name      !<
1810
1811    INTEGER, DIMENSION(1)                             ::  bufshape  !<
1812
1813    REAL(KIND=wp), POINTER, DIMENSION(:)              ::  buf       !<
1814    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:), TARGET ::  data      !<
1815
1816    TYPE(C_PTR)                                       ::  c_data    !<
1817
1818
1819    c_data = C_LOC( data )
1820    bufshape(1) = SIZE( data)
1821    CALL C_F_POINTER( c_data, buf, bufshape )
1822
1823    CALL wrd_mpi_io_global_array_real_1d( name, buf )
1824
1825 END SUBROUTINE wrd_mpi_io_global_array_real_2d
1826
1827
1828
1829!--------------------------------------------------------------------------------------------------!
1830! Description:
1831! ------------
1832!> Write 3d-REAL global array with MPI-IO.
1833!> Array contains identical data on all PEs.
1834!--------------------------------------------------------------------------------------------------!
1835 SUBROUTINE wrd_mpi_io_global_array_real_3d( name, data )
1836
1837    IMPLICIT NONE
1838
1839    CHARACTER(LEN=*), INTENT(IN)                        ::  name      !<
1840
1841    INTEGER, DIMENSION(1)                               ::  bufshape  !<
1842
1843    REAL(KIND=wp), POINTER, DIMENSION(:)                ::  buf       !<
1844    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:), TARGET ::  data      !<
1845
1846    TYPE(C_PTR)                                         ::  c_data    !<
1847
1848
1849    c_data = C_LOC( data )
1850    bufshape(1) = SIZE( data )
1851    CALL C_F_POINTER( c_data, buf, bufshape )
1852
1853    CALL wrd_mpi_io_global_array_real_1d( name, buf )
1854
1855 END SUBROUTINE wrd_mpi_io_global_array_real_3d
1856
1857
1858
1859!--------------------------------------------------------------------------------------------------!
1860! Description:
1861! ------------
1862!> Write 4d-REAL global array with MPI-IO.
1863!> Array contains identical data on all PEs.
1864!--------------------------------------------------------------------------------------------------!
1865 SUBROUTINE wrd_mpi_io_global_array_real_4d( name, data )
1866
1867    IMPLICIT NONE
1868
1869    CHARACTER(LEN=*), INTENT(IN)                          ::  name      !<
1870
1871    INTEGER, DIMENSION(1)                                 ::  bufshape  !<
1872
1873    REAL(KIND=wp), POINTER, DIMENSION(:)                  ::  buf       !<
1874    REAL(KIND=wp), INTENT(IN), DIMENSION(:,:,:,:), TARGET ::  data      !<
1875
1876    TYPE(C_PTR)                                           ::  c_data    !<
1877
1878
1879    c_data = C_LOC( data )
1880    bufshape(1) = SIZE( data)
1881    CALL C_F_POINTER( c_data, buf, bufshape )
1882
1883    CALL wrd_mpi_io_global_array_real_1d( name, buf )
1884
1885 END SUBROUTINE wrd_mpi_io_global_array_real_4d
1886
1887
1888
1889!--------------------------------------------------------------------------------------------------!
1890! Description:
1891! ------------
1892!> Write 1d-INTEGER global array with MPI-IO.
1893!> Array contains identical data on all PEs.
1894!--------------------------------------------------------------------------------------------------!
1895 SUBROUTINE wrd_mpi_io_global_array_int_1d( name, data )
1896
1897    IMPLICIT NONE
1898
1899    CHARACTER(LEN=*), INTENT(IN)                ::  name    !<
1900
1901    INTEGER(KIND=rd_offset_kind)                ::  offset  !<
1902
1903    INTEGER(KIND=iwp), INTENT(IN), DIMENSION(:) ::  data    !<
1904#if defined( __parallel )
1905    INTEGER, DIMENSION(rd_status_size)          ::  status  !<
1906#endif
1907
1908    IF ( header_array_index == max_nr_arrays )  THEN
1909       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
1910    ENDIF
1911
1912    offset = 0
1913    array_names(header_array_index)  = name
1914    array_offset(header_array_index) = array_position
1915    header_array_index = header_array_index + 1
1916
1917!
1918!-- Set default view
1919#if defined( __parallel )
1920    IF ( sm_io%iam_io_pe )  THEN
1921       CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
1922    ENDIF
1923!
1924!-- Only PE 0 writes replicated data
1925    IF ( myid == 0 )  THEN
1926       CALL MPI_FILE_SEEK( fh, array_position, MPI_SEEK_SET, ierr )
1927       CALL MPI_FILE_WRITE( fh, data, SIZE( data), MPI_INTEGER, status, ierr )
1928    ENDIF
1929#else
1930    CALL posix_lseek( fh, array_position )
1931    CALL posix_write( fh, data, SIZE( data ) )
1932#endif
1933    array_position = array_position + SIZE( data ) * 4
1934
1935 END SUBROUTINE wrd_mpi_io_global_array_int_1d
1936
1937
1938
1939!--------------------------------------------------------------------------------------------------!
1940! Description:
1941! ------------
1942!> Read 1d-REAL surface data array with MPI-IO.
1943!--------------------------------------------------------------------------------------------------!
1944 SUBROUTINE rrd_mpi_io_surface( name, data, first_index )
1945
1946    IMPLICIT NONE
1947
1948    CHARACTER(LEN=*), INTENT(IN) ::  name            !<
1949
1950    INTEGER(KIND=rd_offset_kind) ::  disp            !< displacement of actual indices
1951    INTEGER(KIND=rd_offset_kind) ::  disp_f          !< displacement in file
1952    INTEGER(KIND=rd_offset_kind) ::  disp_n          !< displacement of next column
1953    INTEGER(iwp), OPTIONAL       ::  first_index     !<
1954
1955    INTEGER(iwp)                 ::  i               !<
1956    INTEGER(iwp)                 ::  i_f             !<
1957    INTEGER(iwp)                 ::  j               !<
1958    INTEGER(iwp)                 ::  j_f             !<
1959    INTEGER(iwp)                 ::  lo_first_index  !<
1960    INTEGER(iwp)                 ::  nr_bytes        !<
1961    INTEGER(iwp)                 ::  nr_bytes_f      !<
1962    INTEGER(iwp)                 ::  nr_words        !<
1963#if defined( __parallel )
1964    INTEGER, DIMENSION(rd_status_size)  ::  status   !<
1965#endif
1966
1967    LOGICAL                             ::  found    !<
1968
1969    REAL(wp), INTENT(OUT), DIMENSION(:) ::  data     !<
1970
1971
1972    found = .FALSE.
1973    lo_first_index = 1
1974
1975    IF ( MAXVAL( m_global_start ) == -1 )   RETURN   ! Nothing to do on this PE
1976
1977    IF ( PRESENT( first_index ) )   THEN
1978       lo_first_index = first_index
1979    ENDIF
1980
1981    DO  i = 1, tgh%nr_arrays
1982        IF ( TRIM( array_names(i) ) == TRIM( name ) )  THEN
1983           array_position = array_offset(i) + ( lo_first_index - 1 ) *                             &
1984                            total_number_of_surface_values * wp
1985           found = .TRUE.
1986           EXIT
1987        ENDIF
1988    ENDDO
1989
1990    disp   = -1
1991    disp_f = -1
1992    disp_n = -1
1993    IF ( found )  THEN
1994
1995       DO  i = nxl, nxr
1996          DO  j = nys, nyn
1997
1998             IF ( m_global_start(j,i) > 0 )  THEN
1999                disp     = array_position+(m_global_start(j,i)-1) * wp
2000                nr_words = m_end_index(j,i)-m_start_index(j,i)+1
2001                nr_bytes = nr_words * wp
2002             ENDIF
2003             IF ( disp >= 0  .AND.  disp_f == -1 )  THEN   ! First entry
2004                disp_f     = disp
2005                nr_bytes_f = 0
2006                i_f = i
2007                j_f = j
2008             ENDIF
2009             IF ( j == nyn  .AND.  i == nxr )  THEN        ! Last entry
2010                disp_n = -1
2011                IF (  nr_bytes > 0 )  THEN
2012                   nr_bytes_f = nr_bytes_f+nr_bytes
2013                ENDIF
2014             ELSEIF ( j == nyn )  THEN                     ! Next x
2015                IF ( m_global_start(nys,i+1) > 0  .AND.  disp > 0 )  THEN
2016                   disp_n = array_position + ( m_global_start(nys,i+1) - 1 ) * wp
2017                ELSE
2018                   CYCLE
2019                ENDIF
2020             ELSE
2021                IF ( m_global_start(j+1,i) > 0  .AND.  disp > 0 )  THEN
2022                   disp_n = array_position + ( m_global_start(j+1,i) - 1 ) * wp
2023                ELSE
2024                   CYCLE
2025                ENDIF
2026             ENDIF
2027
2028
2029             IF ( disp + nr_bytes == disp_n )  THEN        ! Contiguous block
2030                nr_bytes_f = nr_bytes_f + nr_bytes
2031             ELSE                                          ! Read
2032#if defined( __parallel )
2033                CALL MPI_FILE_SEEK( fhs, disp_f, MPI_SEEK_SET, ierr )
2034                nr_words = nr_bytes_f / wp
2035                CALL MPI_FILE_READ( fhs, data(m_start_index(j_f,i_f)), nr_words, MPI_REAL, status, &
2036                                    ierr )
2037#else
2038                CALL posix_lseek( fh, disp_f )
2039                CALL posix_read( fh, data(m_start_index(j_f,i_f):), nr_bytes_f )
2040#endif
2041                disp_f     = disp
2042                nr_bytes_f = nr_bytes
2043                i_f = i
2044                j_f = j
2045             ENDIF
2046
2047          ENDDO
2048       ENDDO
2049
2050    ELSE
2051
2052       message_string = 'surface array "' // TRIM( name ) // '" not found in restart file'
2053       CALL message( 'rrd_mpi_io_global_array_int_1d', 'PA0722', 3, 2, 0, 6, 0 )
2054
2055    ENDIF
2056
2057!    IF ( lo_first_index == 1 )  THEN
2058!       IF ( debug_level >= 2 .AND. nr_val > 0 )  WRITE(9,*)  'r_surf_1 ', TRIM( name ), ' ', nr_val, SUM( data(1:nr_val) )
2059!    ELSE
2060!       IF ( debug_level >= 2 .AND. nr_val > 0 )  WRITE(9,*)  'r_surf_next ', TRIM( name ), ' ', &
2061!                                                             lo_first_index,nr_val, SUM( data(1:nr_val) )
2062!    ENDIF
2063
2064 END SUBROUTINE rrd_mpi_io_surface
2065
2066
2067
2068!--------------------------------------------------------------------------------------------------!
2069! Description:
2070! ------------
2071!> Read 2d-REAL surface data array with MPI-IO.
2072!> These consist of multiple 1d-REAL surface data arrays.
2073!--------------------------------------------------------------------------------------------------!
2074 SUBROUTINE rrd_mpi_io_surface_2d( name, data )
2075
2076    IMPLICIT NONE
2077
2078    CHARACTER(LEN=*), INTENT(IN)          ::  name  !<
2079
2080    INTEGER(iwp)                          ::  i     !<
2081
2082    REAL(wp), INTENT(OUT), DIMENSION(:,:) ::  data  !<
2083    REAL(wp), DIMENSION(SIZE( data,2))    ::  tmp   !<
2084
2085
2086    DO  i = 1, SIZE( data, 1 )
2087       CALL rrd_mpi_io_surface( name, tmp, first_index = i )
2088       data(i,:) = tmp
2089    ENDDO
2090
2091!
2092!-- Comment from Klaus Ketelsen (September 2018)
2093!-- The intention of the following loop was to let the compiler do the copying on return.
2094!-- In my understanding it is standard conform to pass the second dimension to a 1d-array inside a
2095!-- subroutine and the compiler is responsible to generate code for copying. Acually this works fine
2096!-- for INTENT(IN) variables (wrd_mpi_io_surface_2d). For INTENT(OUT) like in this case the code
2097!-- works on pgi compiler. But both, the Intel 16 and the Cray compiler show wrong answers using
2098!-- this loop. That is the reason why the above auxiliary array tmp was introduced.
2099!    DO  i = 1, SIZE(  data,1)
2100!       CALL rrd_mpi_io_surface( name, data(i,:), first_index = i )
2101!    ENDDO
2102
2103 END SUBROUTINE rrd_mpi_io_surface_2d
2104
2105
2106
2107!--------------------------------------------------------------------------------------------------!
2108! Description:
2109! ------------
2110!> Write 1d-REAL surface data array with MPI-IO.
2111!--------------------------------------------------------------------------------------------------!
2112 SUBROUTINE wrd_mpi_io_surface( name, data, first_index )
2113
2114    IMPLICIT NONE
2115
2116    CHARACTER(LEN=*), INTENT(IN)       ::  name            !<
2117
2118#if defined( __parallel )
2119    INTEGER(KIND=rd_offset_kind)       ::  disp            !<
2120#endif
2121    INTEGER(iwp), OPTIONAL             ::  first_index     !<
2122#if defined( __parallel )
2123    INTEGER(iwp)                       ::  i               !<
2124#endif
2125    INTEGER(iwp)                       ::  lo_first_index  !<
2126    INTEGER(KIND=rd_offset_kind)       ::  offset          !<
2127
2128#if defined( __parallel )
2129    INTEGER, DIMENSION(rd_status_size) ::  status          !<
2130#endif
2131
2132    REAL(wp), INTENT(IN), DIMENSION(:), TARGET ::  data    !<
2133
2134
2135    offset = 0
2136    lo_first_index = 1
2137
2138    IF ( PRESENT( first_index ) )  THEN
2139       lo_first_index = first_index
2140    ENDIF
2141!
2142!-- In case of 2d-data, name is written only once
2143    IF ( lo_first_index == 1 )  THEN
2144
2145       IF ( header_array_index == max_nr_arrays )  THEN
2146          STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
2147       ENDIF
2148
2149       array_names(header_array_index)  = name
2150       array_offset(header_array_index) = array_position
2151       header_array_index = header_array_index + 1
2152
2153    ENDIF
2154
2155#if defined( __parallel )
2156    IF ( sm_io%is_sm_active() )  THEN
2157       DO  i = 1, nr_val
2158          array_1d(i+local_start) = data(i)
2159       ENDDO
2160    ELSE
2161!       array_1d => data                           !kk Did not work in all cases    why???
2162       ALLOCATE( array_1d( SIZE( data ) ) )
2163       array_1d = data
2164    ENDIF
2165
2166    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
2167    IF ( sm_io%iam_io_pe )  THEN
2168       IF ( all_pes_write )  THEN
2169          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_REAL, ft_surf, 'native', MPI_INFO_NULL,  &
2170                                  ierr )
2171          CALL MPI_FILE_WRITE_ALL( fh, array_1d, nr_iope, MPI_REAL, status, ierr )
2172       ELSE
2173          CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2174          IF ( nr_val > 0 )  THEN
2175             disp = array_position + 8 * ( glo_start - 1 )
2176             CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
2177             CALL MPI_FILE_WRITE( fh, array_1d, nr_iope, MPI_REAL, status, ierr )
2178          ENDIF
2179       ENDIF
2180    ENDIF
2181    CALL sm_io%sm_node_barrier()
2182    IF( .NOT. sm_io%is_sm_active() )  DEALLOCATE( array_1d )
2183#else
2184    CALL posix_lseek( fh, array_position )
2185    CALL posix_write( fh, data, nr_val )
2186#endif
2187    array_position = array_position + total_number_of_surface_values * wp
2188
2189!    IF ( lo_first_index == 1 )  THEN
2190!       IF ( debug_level >= 2 .AND. nr_val  > 0 )  WRITE(9,*) 'w_surf_1 ', TRIM( name ), ' ', nr_val, SUM( data(1:nr_val) )
2191!    ELSE
2192!       IF ( debug_level >= 2 .AND. nr_val  > 0 ) WRITE(9,*) 'w_surf_n ', TRIM( name ), ' ', &
2193!                                                            lo_first_index, nr_val, SUM( data(1:nr_val) )
2194!    ENDIF
2195
2196 END SUBROUTINE wrd_mpi_io_surface
2197
2198
2199
2200!--------------------------------------------------------------------------------------------------!
2201! Description:
2202! ------------
2203!> Read 2d-REAL surface data array with MPI-IO.
2204!> This consist of multiple 1d-REAL surface data arrays.
2205!--------------------------------------------------------------------------------------------------!
2206 SUBROUTINE wrd_mpi_io_surface_2d( name, data )
2207
2208    IMPLICIT NONE
2209
2210    CHARACTER(LEN=*), INTENT(IN)         ::  name  !<
2211
2212    INTEGER(iwp)                         ::  i     !<
2213
2214    REAL(wp), INTENT(IN), DIMENSION(:,:) ::  data  !<
2215
2216
2217    DO  i = 1, SIZE( data, 1 )
2218       CALL wrd_mpi_io_surface( name, data(i,:), first_index = i )
2219    ENDDO
2220
2221 END SUBROUTINE wrd_mpi_io_surface_2d
2222
2223
2224
2225!--------------------------------------------------------------------------------------------------!
2226! Description:
2227! ------------
2228!> Close restart file for MPI-IO
2229!--------------------------------------------------------------------------------------------------!
2230 SUBROUTINE rd_mpi_io_close
2231
2232    IMPLICIT NONE
2233
2234    INTEGER(iwp)                       ::  gh_size  !<
2235    INTEGER(KIND=rd_offset_kind)       ::  offset   !<
2236#if defined( __parallel )
2237    INTEGER, DIMENSION(rd_status_size) ::  status   !<
2238#endif
2239
2240#if ! defined( __parallel )
2241    TYPE(C_PTR)                        ::  buf_ptr  !<
2242#endif
2243
2244
2245    offset = 0
2246
2247    IF ( wr_flag  .AND.  sm_io%iam_io_pe )  THEN
2248
2249       tgh%nr_int    = header_int_index - 1
2250       tgh%nr_char   = header_char_index - 1
2251       tgh%nr_real   = header_real_index - 1
2252       tgh%nr_arrays = header_array_index - 1
2253       tgh%total_nx  = lb%nx + 1
2254       tgh%total_ny  = lb%ny + 1
2255       IF ( include_total_domain_boundaries )  THEN   ! Not sure, if LOGICAL interpretation is the same for all compilers,
2256          tgh%i_outer_bound = 1                       ! therefore store as INTEGER in general header
2257       ELSE
2258          tgh%i_outer_bound = 0
2259       ENDIF
2260!
2261!--    Check for big/little endian format. This check is currently not used, and could be removed
2262!--    if we can assume little endian as the default on all machines.
2263       CALL rd_mpi_io_check_endian( tgh%endian )
2264
2265!
2266!--    Set default view
2267#if defined( __parallel )
2268       CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2269#endif
2270!
2271!--    Write header file
2272       gh_size = storage_size(tgh) / 8
2273       IF ( myid == 0 )  THEN   ! myid = 0 always performs I/O, even if I/O is limited to some cores
2274#if defined( __parallel )
2275          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2276          CALL MPI_FILE_WRITE( fh, tgh, gh_size, MPI_INT, status, ierr )
2277          header_position = header_position + gh_size
2278!
2279!--       INTEGER values
2280          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2281          CALL MPI_FILE_WRITE( fh, int_names, SIZE( int_names ) * 32, MPI_CHAR, status, ierr )
2282          header_position = header_position + SIZE( int_names ) * 32
2283
2284          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2285          CALL MPI_FILE_WRITE( fh, int_values, SIZE( int_values ), MPI_INT, status, ierr )
2286          header_position = header_position + SIZE( int_values ) * iwp
2287!
2288!--       Character entries
2289          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2290          CALL MPI_FILE_WRITE( fh, text_lines, SIZE( text_lines ) * 128, MPI_CHAR, status, ierr )
2291          header_position = header_position + SIZE( text_lines ) * 128
2292!
2293!---      REAL values
2294          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2295          CALL MPI_FILE_WRITE( fh, real_names, SIZE( real_names ) * 32, MPI_CHAR, status, ierr )
2296          header_position = header_position + SIZE( real_names ) * 32
2297
2298          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2299          CALL MPI_FILE_WRITE( fh, real_values, SIZE( real_values ), MPI_REAL, status, ierr )
2300          header_position = header_position + SIZE( real_values ) * wp
2301!
2302!--       2d- and 3d- distributed array headers, all replicated array headers
2303          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2304          CALL MPI_FILE_WRITE( fh, array_names, SIZE( array_names ) * 32, MPI_CHAR, status, ierr )
2305          header_position = header_position + SIZE( array_names ) * 32
2306
2307          CALL MPI_FILE_SEEK( fh, header_position, MPI_SEEK_SET, ierr )
2308          CALL MPI_FILE_WRITE( fh, array_offset, SIZE( array_offset ) * MPI_OFFSET_KIND, MPI_BYTE, &
2309                               status, ierr )  ! There is no I*8 datatype in Fortran
2310          header_position = header_position + SIZE( array_offset ) * rd_offset_kind
2311#else
2312          CALL posix_lseek( fh, header_position )
2313          buf_ptr = C_LOC( tgh )
2314          CALL posix_write( fh, buf_ptr, gh_size )
2315          header_position = header_position + gh_size
2316!
2317!--       INTEGER values
2318          CALL posix_lseek( fh, header_position )
2319          CALL posix_write( fh, int_names )
2320          header_position = header_position + SIZE( int_names ) * 32
2321
2322          CALL posix_lseek( fh, header_position )
2323          CALL posix_write( fh, int_values, SIZE( int_values ) )
2324          header_position = header_position + SIZE( int_values ) * iwp
2325!
2326!--       Character entries
2327          CALL posix_lseek( fh, header_position )
2328          CALL posix_write( fh, text_lines )
2329          header_position = header_position + SIZE( text_lines ) * 128
2330!
2331!--       REAL values
2332          CALL posix_lseek( fh, header_position )
2333          CALL posix_write( fh, real_names )
2334          header_position = header_position + SIZE( real_names ) * 32
2335
2336          CALL posix_lseek( fh, header_position )
2337          CALL posix_write( fh, real_values, SIZE( real_values ) )
2338          header_position = header_position + SIZE( real_values ) * wp
2339!
2340!--       2d- and 3d-distributed array headers, all replicated array headers
2341          CALL posix_lseek( fh, header_position )
2342          CALL posix_write( fh, array_names )
2343          header_position = header_position + SIZE( array_names ) * 32
2344
2345          CALL posix_lseek( fh, header_position )
2346          CALL posix_write( fh, array_offset, SIZE( array_offset ) )
2347          header_position = header_position + SIZE( array_offset ) * rd_offset_kind
2348#endif
2349          IF ( debug_output )  CALL rd_mpi_io_print_header
2350       ENDIF
2351
2352    ENDIF
2353
2354!
2355!-- Free file types
2356    CALL rd_mpi_io_free_filetypes
2357
2358!
2359!-- Close MPI-IO files
2360#if defined( __parallel )
2361!
2362!-- Restart file has been opened with comm2d
2363    IF ( fhs /= -1 )  THEN
2364       CALL MPI_FILE_CLOSE( fhs, ierr )
2365    ENDIF
2366#endif
2367
2368    IF ( sm_io%iam_io_pe )  THEN
2369
2370#if defined( __parallel )
2371       CALL MPI_FILE_CLOSE( fh, ierr )
2372#else
2373       CALL posix_close( fh )
2374#endif
2375
2376    ENDIF
2377
2378    restart_file_size = array_position / ( 1024.0_dp * 1024.0_dp )
2379
2380 END SUBROUTINE rd_mpi_io_close
2381
2382
2383
2384!--------------------------------------------------------------------------------------------------!
2385! Description:
2386! ------------
2387!> This subroutine prepares a filetype and some variables for the I/O of surface data.
2388!> Whenever a new set of start_index and end_index is used, rd_mpi_io_surface_filetypes has to be
2389!> called. A main feature of this subroutine is computing the global start indices of the 1d- and
2390!> 2d- surface arrays.
2391!> Even if I/O is done by a limited number of cores only, the surface data are read by ALL cores!
2392!> Reading them by some cores and then distributing the data would result in complicated code
2393!> which is suspectable for errors and overloads the reading subroutine. Since reading of surface
2394!> data is not time critical (data size is comparably small), it will be read by all cores.
2395!--------------------------------------------------------------------------------------------------!
2396 SUBROUTINE rd_mpi_io_surface_filetypes( start_index, end_index, data_to_write, global_start )
2397
2398    IMPLICIT NONE
2399
2400    INTEGER(iwp)                          ::  i           !<  loop index
2401    INTEGER(KIND=rd_offset_kind)          ::  offset      !<
2402
2403    INTEGER(iwp), DIMENSION(1)            ::  dims1       !<
2404    INTEGER(iwp), DIMENSION(1)            ::  lize1       !<
2405    INTEGER(iwp), DIMENSION(1)            ::  start1      !<
2406
2407    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  all_nr_val  !< number of values for all PEs
2408    INTEGER(iwp), DIMENSION(0:numprocs-1) ::  lo_nr_val   !< local number of values in x and y direction
2409
2410
2411    INTEGER, INTENT(IN), DIMENSION(nys:nyn,nxl:nxr)  ::  end_index     !<
2412    INTEGER, INTENT(OUT), DIMENSION(nys:nyn,nxl:nxr) ::  global_start  !<
2413    INTEGER, INTENT(IN), DIMENSION(nys:nyn,nxl:nxr)  ::  start_index   !<
2414
2415    LOGICAL, INTENT(OUT) ::  data_to_write  !< returns, if surface data have to be written
2416
2417
2418    offset = 0
2419    lo_nr_val= 0
2420    lo_nr_val(myid) = MAXVAL( end_index )
2421#if defined( __parallel )
2422    CALL MPI_ALLREDUCE( lo_nr_val, all_nr_val, numprocs, MPI_INTEGER, MPI_SUM, comm2d, ierr )
2423    IF ( ft_surf /= -1  .AND.  sm_io%iam_io_pe )  THEN
2424       CALL MPI_TYPE_FREE( ft_surf, ierr )    ! If set, free last surface filetype
2425    ENDIF
2426
2427    IF ( win_surf /= -1 )  THEN
2428       IF ( sm_io%is_sm_active() )  THEN
2429          CALL MPI_WIN_FREE( win_surf, ierr )
2430       ENDIF
2431       win_surf = -1
2432    ENDIF
2433
2434    IF ( sm_io%is_sm_active() .AND. rd_flag )  THEN
2435       IF ( fhs == -1 )  THEN
2436          CALL MPI_FILE_OPEN( comm2d, TRIM( io_file_name ), MPI_MODE_RDONLY, MPI_INFO_NULL, fhs,   &
2437                              ierr )
2438       ENDIF
2439    ELSE
2440       fhs = fh
2441    ENDIF
2442#else
2443    all_nr_val(myid) = lo_nr_val(myid)
2444#endif
2445    nr_val = lo_nr_val(myid)
2446
2447    total_number_of_surface_values = 0
2448    DO  i = 0, numprocs-1
2449       IF ( i == myid )  THEN
2450          glo_start = total_number_of_surface_values + 1
2451       ENDIF
2452       total_number_of_surface_values = total_number_of_surface_values + all_nr_val(i)
2453    ENDDO
2454
2455!
2456!-- Actions during reading
2457    IF ( rd_flag )  THEN
2458       IF ( .NOT. ALLOCATED( m_start_index )  )  ALLOCATE( m_start_index(nys:nyn,nxl:nxr)  )
2459       IF ( .NOT. ALLOCATED( m_end_index )    )  ALLOCATE( m_end_index(nys:nyn,nxl:nxr)    )
2460       IF ( .NOT. ALLOCATED( m_global_start ) )  ALLOCATE( m_global_start(nys:nyn,nxl:nxr) )
2461!
2462!--    Save arrays for later reading
2463       m_start_index  = start_index
2464       m_end_index    = end_index
2465       m_global_start = global_start
2466       nr_val = MAXVAL( end_index )
2467
2468#if defined( __parallel )
2469       CALL MPI_FILE_SET_VIEW( fhs, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
2470#endif
2471    ENDIF
2472
2473!
2474!-- Actions during writing
2475    IF ( wr_flag )  THEN
2476!
2477!--    Create surface filetype
2478       ft_surf      = -1
2479       global_start = start_index + glo_start - 1
2480
2481       WHERE ( end_index < start_index )
2482          global_start = -1
2483       ENDWHERE
2484
2485#if defined( __parallel )
2486       IF ( sm_io%is_sm_active() )  THEN
2487          IF ( sm_io%iam_io_pe )  THEN
2488!
2489!--          Calculate number of values of all PEs of an I/O group
2490             nr_iope = 0
2491             DO  i = myid, myid+sm_io%sh_npes-1
2492                nr_iope = nr_iope + all_nr_val(i)
2493             ENDDO
2494          ELSE
2495             local_start = 0
2496             DO  i = myid-sm_io%sh_rank, myid-1
2497                local_start = local_start + all_nr_val(i)
2498             ENDDO
2499          ENDIF
2500!
2501!--       Get the size of shared memory window on all PEs
2502          CALL MPI_BCAST( nr_iope, 1, MPI_INTEGER, 0, sm_io%comm_shared, ierr )
2503          CALL sm_io%sm_allocate_shared( array_1d, 1, MAX( 1, nr_iope ), win_surf )
2504       ELSE
2505          nr_iope = nr_val
2506       ENDIF
2507#else
2508       nr_iope = nr_val
2509#endif
2510
2511!
2512!--    Check, if surface data exist on this PE
2513       data_to_write = .TRUE.
2514       IF ( total_number_of_surface_values == 0 )  THEN
2515          data_to_write = .FALSE.
2516          RETURN
2517       ENDIF
2518
2519       IF ( sm_io%iam_io_pe )  THEN
2520
2521          all_pes_write = ( MINVAL( all_nr_val ) > 0 )
2522
2523          IF ( all_pes_write )  THEN
2524             dims1(1)  = total_number_of_surface_values
2525             lize1(1)  = nr_iope
2526             start1(1) = glo_start-1
2527
2528#if defined( __parallel )
2529             IF ( total_number_of_surface_values > 0 )  THEN
2530                 CALL MPI_TYPE_CREATE_SUBARRAY( 1, dims1, lize1, start1, MPI_ORDER_FORTRAN,        &
2531                                                MPI_REAL, ft_surf, ierr )
2532                 CALL MPI_TYPE_COMMIT( ft_surf, ierr )
2533             ENDIF
2534#endif
2535          ENDIF
2536       ENDIF
2537
2538    ENDIF
2539
2540 END SUBROUTINE rd_mpi_io_surface_filetypes
2541
2542
2543
2544!--------------------------------------------------------------------------------------------------!
2545! Description:
2546! ------------
2547!> This subroutine creates file types to access 2d-/3d-REAL arrays and 2d-INTEGER arrays
2548!> distributed in blocks among processes to a single file that contains the global arrays.
2549!--------------------------------------------------------------------------------------------------!
2550  SUBROUTINE rd_mpi_io_create_filetypes
2551
2552    IMPLICIT NONE
2553
2554    INTEGER, DIMENSION(2) ::  dims2   !<
2555    INTEGER, DIMENSION(2) ::  lize2   !<
2556    INTEGER, DIMENSION(2) ::  start2  !<
2557
2558    INTEGER, DIMENSION(3) ::  dims3   !<
2559    INTEGER, DIMENSION(3) ::  lize3   !<
2560    INTEGER, DIMENSION(3) ::  start3  !<
2561
2562    TYPE(local_boundaries) ::  save_io_grid  !< temporary variable to store grid settings
2563
2564
2565    IF ( sm_io%is_sm_active() )  THEN
2566       save_io_grid = sm_io%io_grid
2567    ENDIF
2568
2569!
2570!-- Decision, if storage with or without ghost layers.
2571!-- Please note that the indexing of the global array always starts at 0, even in Fortran.
2572!-- Therefore the PE local indices have to be shifted by nbgp in the case with ghost layers.
2573    IF ( include_total_domain_boundaries )  THEN
2574
2575       lb%nxl = nxl + nbgp
2576       lb%nxr = nxr + nbgp
2577       lb%nys = nys + nbgp
2578       lb%nyn = nyn + nbgp
2579       lb%nnx = nnx
2580       lb%nny = nny
2581       lb%nx  = nx + 2 * nbgp
2582       lb%ny  = ny + 2 * nbgp
2583       IF ( myidx == 0 )  THEN
2584          lb%nxl = lb%nxl - nbgp
2585          lb%nnx = lb%nnx + nbgp
2586       ENDIF
2587       IF ( myidx == npex-1  .OR.  npex == -1 )  THEN   ! npex == 1 if -D__parallel not set
2588          lb%nxr = lb%nxr + nbgp
2589          lb%nnx = lb%nnx + nbgp
2590       ENDIF
2591       IF ( myidy == 0 )  THEN
2592          lb%nys = lb%nys - nbgp
2593          lb%nny = lb%nny + nbgp
2594       ENDIF
2595       IF ( myidy == npey-1  .OR.  npey == -1 )  THEN   ! npey == 1 if -D__parallel not set
2596          lb%nyn = lb%nyn + nbgp
2597          lb%nny = lb%nny + nbgp
2598       ENDIF
2599
2600       CALL sm_io%sm_adjust_outer_boundary()
2601
2602    ELSE
2603
2604       lb%nxl = nxl
2605       lb%nxr = nxr
2606       lb%nys = nys
2607       lb%nyn = nyn
2608       lb%nnx = nnx
2609       lb%nny = nny
2610       lb%nx  = nx
2611       lb%ny  = ny
2612
2613    ENDIF
2614
2615    IF ( sm_io%is_sm_active() )  THEN
2616#if defined( __parallel )
2617       CALL sm_io%sm_allocate_shared( array_2d,  sm_io%io_grid%nxl, sm_io%io_grid%nxr,             &
2618                                      sm_io%io_grid%nys, sm_io%io_grid%nyn, win_2dr )
2619       CALL sm_io%sm_allocate_shared( array_2di, save_io_grid%nxl, save_io_grid%nxr,               &
2620                                      save_io_grid%nys, save_io_grid%nyn, win_2di )
2621       CALL sm_io%sm_allocate_shared( array_3d, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,  &
2622                                      sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3dr )
2623#endif
2624    ELSE
2625       ALLOCATE( array_2d(lb%nxl:lb%nxr,lb%nys:lb%nyn) )
2626       ALLOCATE( array_2di(nxl:nxr,nys:nyn) )
2627       ALLOCATE( array_3d(nzb:nzt+1,lb%nxl:lb%nxr,lb%nys:lb%nyn) )
2628       sm_io%io_grid = lb
2629    ENDIF
2630
2631!
2632!-- Create filetype for 2d-REAL array with ghost layers around the total domain
2633    dims2(1)  = lb%nx + 1
2634    dims2(2)  = lb%ny + 1
2635
2636    lize2(1)  = sm_io%io_grid%nnx
2637    lize2(2)  = sm_io%io_grid%nny
2638
2639    start2(1) = sm_io%io_grid%nxl
2640    start2(2) = sm_io%io_grid%nys
2641
2642#if defined( __parallel )
2643    IF ( sm_io%iam_io_pe )  THEN
2644       CALL MPI_TYPE_CREATE_SUBARRAY( 2, dims2, lize2, start2, MPI_ORDER_FORTRAN, MPI_REAL,        &
2645                                      ft_2d, ierr )
2646       CALL MPI_TYPE_COMMIT( ft_2d, ierr )
2647    ENDIF
2648#endif
2649!
2650!-- Create filetype for 2d-INTEGER array without ghost layers around the total domain
2651    dims2(1)  = nx + 1
2652    dims2(2)  = ny + 1
2653
2654    IF ( sm_io%is_sm_active() )  THEN
2655
2656       lize2(1)  = save_io_grid%nnx
2657       lize2(2)  = save_io_grid%nny
2658
2659       start2(1) = save_io_grid%nxl
2660       start2(2) = save_io_grid%nys
2661
2662    ELSE
2663
2664       lize2(1)  = nnx
2665       lize2(2)  = nny
2666
2667       start2(1) = nxl
2668       start2(2) = nys
2669
2670    ENDIF
2671
2672#if defined( __parallel )
2673    IF ( sm_io%iam_io_pe )  THEN
2674       CALL MPI_TYPE_CREATE_SUBARRAY( 2, dims2, lize2, start2, MPI_ORDER_FORTRAN, MPI_INTEGER,     &
2675                                      ft_2di_nb, ierr )
2676       CALL MPI_TYPE_COMMIT( ft_2di_nb, ierr )
2677    ENDIF
2678#endif
2679!
2680!-- Create filetype for 3d-REAL array
2681    dims3(1)  = nz + 2
2682    dims3(2)  = lb%nx + 1
2683    dims3(3)  = lb%ny + 1
2684
2685    lize3(1)  = dims3(1)
2686    lize3(2)  = sm_io%io_grid%nnx
2687    lize3(3)  = sm_io%io_grid%nny
2688
2689    start3(1) = nzb
2690    start3(2) = sm_io%io_grid%nxl
2691    start3(3) = sm_io%io_grid%nys
2692
2693#if defined( __parallel )
2694    IF ( sm_io%iam_io_pe )  THEN
2695       CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_REAL, ft_3d, &
2696                                      ierr )
2697       CALL MPI_TYPE_COMMIT( ft_3d, ierr )
2698    ENDIF
2699#endif
2700
2701 END SUBROUTINE rd_mpi_io_create_filetypes
2702
2703
2704
2705!--------------------------------------------------------------------------------------------------!
2706! Description:
2707! ------------
2708!> This subroutine creates file types to access 3d-soil arrays distributed in blocks among processes
2709!> to a single file that contains the global arrays. It is not required for the serial mode.
2710!--------------------------------------------------------------------------------------------------!
2711#if defined( __parallel )
2712 SUBROUTINE rd_mpi_io_create_filetypes_3dsoil( nzb_soil, nzt_soil )
2713
2714    IMPLICIT NONE
2715
2716    INTEGER, INTENT(IN)   ::  nzb_soil  !<
2717    INTEGER, INTENT(IN)   ::  nzt_soil  !<
2718
2719    INTEGER, DIMENSION(3) ::  dims3     !<
2720    INTEGER, DIMENSION(3) ::  lize3     !<
2721    INTEGER, DIMENSION(3) ::  start3    !<
2722
2723
2724    IF ( sm_io%is_sm_active() )  THEN
2725       CALL sm_io%sm_allocate_shared( array_3d_soil, nzb_soil, nzt_soil, sm_io%io_grid%nxl,        &
2726                                      sm_io%io_grid%nxr, sm_io%io_grid%nys, sm_io%io_grid%nyn,     &
2727                                      win_3ds )
2728    ELSE
2729       ALLOCATE( array_3d_soil(nzb_soil:nzt_soil,lb%nxl:lb%nxr,lb%nys:lb%nyn) )
2730       sm_io%io_grid = lb
2731    ENDIF
2732
2733!
2734!-- Create filetype for 3d-soil array
2735    dims3(1)  = nzt_soil - nzb_soil + 1
2736    dims3(2)  = lb%nx + 1
2737    dims3(3)  = lb%ny + 1
2738
2739    lize3(1)  = dims3(1)
2740    lize3(2)  = sm_io%io_grid%nnx
2741    lize3(3)  = sm_io%io_grid%nny
2742
2743    start3(1) = nzb_soil
2744    start3(2) = sm_io%io_grid%nxl
2745    start3(3) = sm_io%io_grid%nys
2746
2747    IF ( sm_io%iam_io_pe )  THEN
2748       CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_REAL,        &
2749                                      ft_3dsoil, ierr )
2750       CALL MPI_TYPE_COMMIT( ft_3dsoil, ierr )
2751    ENDIF
2752
2753 END SUBROUTINE rd_mpi_io_create_filetypes_3dsoil
2754#endif
2755
2756
2757
2758!--------------------------------------------------------------------------------------------------!
2759! Description:
2760! ------------
2761!> Free all file types that have been created for MPI-IO.
2762!--------------------------------------------------------------------------------------------------!
2763 SUBROUTINE rd_mpi_io_free_filetypes
2764
2765    IMPLICIT NONE
2766
2767
2768#if defined( __parallel )
2769    IF ( filetypes_created )  THEN
2770
2771       IF ( sm_io%iam_io_pe )  THEN
2772          CALL MPI_TYPE_FREE( ft_2d, ierr )
2773          CALL MPI_TYPE_FREE( ft_2di_nb, ierr )
2774          CALL MPI_TYPE_FREE( ft_3d, ierr )
2775       ENDIF
2776
2777       IF ( sm_io%is_sm_active() )  THEN
2778          CALL sm_io%sm_free_shared( win_2dr )
2779          CALL sm_io%sm_free_shared( win_2di )
2780          CALL sm_io%sm_free_shared( win_3dr )
2781       ELSE
2782          DEALLOCATE( array_2d, array_2di, array_3d )
2783       ENDIF
2784
2785    ENDIF
2786!
2787!-- Free last surface filetype
2788    IF ( sm_io%iam_io_pe .AND. ft_surf /= -1 )  THEN
2789       CALL MPI_TYPE_FREE( ft_surf, ierr )
2790    ENDIF
2791
2792    IF ( sm_io%is_sm_active() .AND.  win_surf /= -1 )  THEN
2793       CALL sm_io%sm_free_shared( win_surf )
2794    ENDIF
2795
2796    ft_surf  = -1
2797    win_surf = -1
2798#else
2799    DEALLOCATE( array_2d, array_2di, array_3d )
2800#endif
2801
2802 END SUBROUTINE rd_mpi_io_free_filetypes
2803
2804
2805
2806!--------------------------------------------------------------------------------------------------!
2807! Description:
2808! ------------
2809!> Print the restart data file header (MPI-IO format) for debugging.
2810!--------------------------------------------------------------------------------------------------!
2811 SUBROUTINE rd_mpi_io_print_header
2812
2813    IMPLICIT NONE
2814
2815    INTEGER(iwp) ::  i  !<
2816
2817
2818    WRITE (9,*)  'header position after reading the restart file header: ', header_position
2819    WRITE (9,*)  ' '
2820    WRITE (9,*)  'restart file header content:'
2821    WRITE (9,*)  '----------------------------'
2822    WRITE (9,*)  ' '
2823
2824    WRITE (9,*)  ' CHARACTER header values   Total number: ', tgh%nr_char
2825    WRITE (9,*)  ' '
2826    DO  i = 1, tgh%nr_char
2827       WRITE( 9, '(I3,A,1X,A)' )  i, ': ', text_lines(i)(1:80)
2828    ENDDO
2829    WRITE (9,*)  ' '
2830
2831    WRITE (9,*) ' INTEGER header variables and values   Total number: ', tgh%nr_int
2832    WRITE (9,*)  ' '
2833    DO  i = 1, tgh%nr_int
2834       WRITE(9,*)  ' variable: ', int_names(i), '  value: ', int_values(i)
2835    ENDDO
2836    WRITE (9,*)  ' '
2837
2838    WRITE (9,*)  ' REAL header variables and values   Total number: ', tgh%nr_real
2839    WRITE (9,*)  ' '
2840    DO  i = 1, tgh%nr_real
2841       WRITE(9,*)  ' variable: ', real_names(i), '  value: ', real_values(i)
2842    ENDDO
2843    WRITE (9,*)  ' '
2844
2845    WRITE (9,*)  ' Header entries with offset (2d/3d arrays)   Total number: ', tgh%nr_arrays
2846    WRITE (9,*)  ' '
2847    DO  i = 1, tgh%nr_arrays
2848       WRITE(9,*)  ' variable: ', array_names(i), '  offset: ', array_offset(i)
2849    ENDDO
2850    WRITE (9,*)  ' '
2851
2852 END SUBROUTINE rd_mpi_io_print_header
2853
2854
2855
2856!--------------------------------------------------------------------------------------------------!
2857! Description:
2858! ------------
2859!> Check if big/little endian data format is used.
2860!> An int*4 pointer is set to a int*8 variable, the int*8 is set to 1, and then it is checked, if
2861!> the first 4 bytes of the pointer are equal 1 (little endian) or not.
2862!--------------------------------------------------------------------------------------------------!
2863 SUBROUTINE rd_mpi_io_check_endian( i_endian )
2864
2865    IMPLICIT NONE
2866
2867    INTEGER, INTENT(out)                   ::  i_endian  !<
2868    INTEGER(KIND=8), TARGET                ::  int8      !<
2869
2870    INTEGER, DIMENSION(1)                  ::  bufshape  !<
2871    INTEGER(KIND=4), POINTER, DIMENSION(:) ::  int4      !<
2872
2873    TYPE(C_PTR)                            ::  ptr       !<
2874
2875
2876    ptr = C_LOC( int8 )
2877    bufshape(1) = 2
2878    CALL C_F_POINTER( ptr, int4, bufshape )
2879
2880    int8 = 1
2881
2882    IF ( int4(1) == 1 )  THEN
2883       i_endian = 1    ! Little endian
2884    ELSE
2885       i_endian = 2    ! Big endian
2886    ENDIF
2887
2888 END SUBROUTINE rd_mpi_io_check_endian
2889
2890 END MODULE restart_data_mpi_io_mod
Note: See TracBrowser for help on using the repository browser.