Home


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

Ignore:
Timestamp:
Jul 29, 2020 7:23:03 AM (4 years ago)
Author:
raasch
Message:

extensions required for MPI-I/O of particle data to restart files

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/restart_data_mpi_io_mod.f90

    r4622 r4628  
    2525! -----------------
    2626! $Id$
     27! extensions required for MPI-I/O of particle data to restart files
     28!
     29! 4622 2020-07-23 09:02:23Z raasch
    2730! unused variables removed
    2831!
    … …  
    6568!
    6669! 4495 2020-04-13 20:11:20Z raasch
    67 ! Initial version (K. Ketelsen), adjusted to PALM formatting standards (s. Raasch)
    68 !
     70! Initial version (K. Ketelsen), adjusted to PALM formatting standards (S. Raasch)
    6971!
    7072!
    … …  
    127129    USE kinds
    128130
     131    USE particle_attributes,                                                                       &
     132        ONLY:  grid_particles,                                                                     &
     133               particles,                                                                          &
     134               particle_type,                                                                      &
     135               prt_count,                                                                          &
     136               zero_particle
     137
    129138    USE pegrid,                                                                                    &
    130139        ONLY:  comm1dx,                                                                            &
    … …  
    161170
    162171    INTEGER(iwp)            ::  comm_io          !< Communicator for MPI-IO
    163     INTEGER(iwp)            ::  fh               !< MPI-IO file handle
     172    INTEGER(iwp)            ::  fh = -1          !< MPI-IO file handle
    164173#if defined( __parallel )
    165174    INTEGER(iwp)            ::  fhs = -1         !< MPI-IO file handle to open file with comm2d always
    … …  
    170179    INTEGER(iwp)            ::  ft_2d            !< MPI filetype 2D array REAL with outer boundaries
    171180    INTEGER(iwp)            ::  ft_3d            !< MPI filetype 3D array REAL with outer boundaries
     181    INTEGER(iwp)            ::  ft_3di4 = -1     !< MPI filetype 3D array INTEGER*4
     182    INTEGER(iwp)            ::  ft_3di8 = -1     !< MPI filetype 3D array INTEGER*8
    172183    INTEGER(iwp)            ::  ft_3dsoil        !< MPI filetype for 3d-soil array
    173184#endif
    … …  
    181192    INTEGER(iwp)            ::  win_2di          !<
    182193    INTEGER(iwp)            ::  win_2dr          !<
     194    INTEGER(iwp)            ::  win_3di4 = -1    !<
     195    INTEGER(iwp)            ::  win_3di8 = -1    !<
    183196    INTEGER(iwp)            ::  win_3dr          !<
    184197    INTEGER(iwp)            ::  win_3ds          !<
    … …  
    190203    INTEGER(KIND=rd_offset_kind) ::  header_position  !<
    191204
    192     INTEGER(iwp), DIMENSION(:,:), POINTER, CONTIGUOUS ::  array_2di  !<
     205    INTEGER(iwp), DIMENSION(:,:), POINTER, CONTIGUOUS   ::  array_2di   !<
    193206
    194207    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_end_index     !<
    … …  
    196209    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  m_start_index   !<
    197210
     211    INTEGER(isp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  array_3di4  !<
     212    INTEGER(idp), DIMENSION(:,:,:), POINTER, CONTIGUOUS ::  array_3di8  !<
    198213
    199214    LOGICAL ::  all_pes_write                 !< all PEs have data to write
    … …  
    279294    TYPE(domain_decomposition_grid_features) ::  prerun_grid   !< grid variables for the prerun
    280295
     296!
     297!-- MPI_INTEGER8 is not standard MPI, but is supported on most MPI distibutions
     298!-- If not suppported, a workaround could be enabled with the following preprocessor directive
     299!#if defined( __NO_INTEGER8)
     300!    INTEGER ::  MPI_INTEGER8  !< MPI data type INTEGER8
     301!#endif
    281302
    282303    SAVE
    … …  
    296317    END INTERFACE rd_mpi_io_close
    297318
     319    INTERFACE rd_mpi_io_check_open
     320       MODULE PROCEDURE rd_mpi_io_check_open
     321    END INTERFACE rd_mpi_io_check_open
     322
    298323    INTERFACE rd_mpi_io_open
    299324       MODULE PROCEDURE rd_mpi_io_open
    … …  
    304329       MODULE PROCEDURE rrd_mpi_io_int
    305330       MODULE PROCEDURE rrd_mpi_io_int_2d
     331       MODULE PROCEDURE rrd_mpi_io_int4_3d
     332       MODULE PROCEDURE rrd_mpi_io_int8_3d
    306333       MODULE PROCEDURE rrd_mpi_io_logical
    307334       MODULE PROCEDURE rrd_mpi_io_real
    … …  
    324351    END INTERFACE rrd_mpi_io_surface
    325352
     353    INTERFACE rrd_mpi_io_particles
     354       MODULE PROCEDURE rrd_mpi_io_particles
     355    END INTERFACE rrd_mpi_io_particles
     356
     357    INTERFACE rd_mpi_io_particle_filetypes
     358       MODULE PROCEDURE rd_mpi_io_particle_filetypes
     359    END INTERFACE rd_mpi_io_particle_filetypes
     360
    326361    INTERFACE rd_mpi_io_surface_filetypes
    327362       MODULE PROCEDURE rd_mpi_io_surface_filetypes
    … …  
    332367       MODULE PROCEDURE wrd_mpi_io_int
    333368       MODULE PROCEDURE wrd_mpi_io_int_2d
     369       MODULE PROCEDURE wrd_mpi_io_int4_3d
     370       MODULE PROCEDURE wrd_mpi_io_int8_3d
    334371       MODULE PROCEDURE wrd_mpi_io_logical
    335372       MODULE PROCEDURE wrd_mpi_io_real
    … …  
    347384    END INTERFACE wrd_mpi_io_global_array
    348385
     386    INTERFACE wrd_mpi_io_particles
     387       MODULE PROCEDURE wrd_mpi_io_particles
     388    END INTERFACE wrd_mpi_io_particles
     389
    349390    INTERFACE wrd_mpi_io_surface
    350391       MODULE PROCEDURE wrd_mpi_io_surface
    … …  
    353394
    354395    PUBLIC  rd_mpi_io_check_array,                                                                 &
     396            rd_mpi_io_check_open,                                                                  &
    355397            rd_mpi_io_close,                                                                       &
    356398            rd_mpi_io_open,                                                                        &
     399            rd_mpi_io_particle_filetypes,                                                          &
     400            rd_mpi_io_surface_filetypes,                                                           &
    357401            rrd_mpi_io,                                                                            &
    358402            rrd_mpi_io_global_array,                                                               &
     403            rrd_mpi_io_particles,                                                                  &
    359404            rrd_mpi_io_surface,                                                                    &
    360             rd_mpi_io_surface_filetypes,                                                           &
    361405            wrd_mpi_io,                                                                            &
    362406            wrd_mpi_io_global_array,                                                               &
     407            wrd_mpi_io_particles,                                                                  &
    363408            wrd_mpi_io_surface
    364409
    … …  
    394439    TYPE(C_PTR)                   ::  buf_ptr  !<
    395440#endif
    396 
    397 !    write(9,*) 'Here is rd_mpi_io_open',nx,nx_on_file,ny,ny_on_file,TRIM(action)   !kk may become Debug Output
    398441
    399442    offset = 0
    … …  
    14151458! Description:
    14161459! ------------
     1460!> Read 3d-INTEGER*4 array with MPI-IO
     1461!--------------------------------------------------------------------------------------------------!
     1462 SUBROUTINE rrd_mpi_io_int4_3d( name, data )
     1463
     1464    IMPLICIT NONE
     1465
     1466    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     1467
     1468    INTEGER(iwp)                       ::  i       !<
     1469
     1470#if defined( __parallel )
     1471    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     1472#endif
     1473
     1474    LOGICAL                            ::  found   !<
     1475
     1476    INTEGER(isp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data  !<
     1477
     1478
     1479    found = .FALSE.
     1480    data  = -1.0
     1481
     1482    DO  i = 1, tgh%nr_arrays
     1483       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
     1484          array_position = array_offset(i)
     1485          found = .TRUE.
     1486          EXIT
     1487       ENDIF
     1488    ENDDO
     1489
     1490    IF ( found )  THEN
     1491
     1492#if defined( __parallel )
     1493       CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited # of cores is inactive
     1494       IF( sm_io%iam_io_pe )  THEN
     1495          CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_3di4, 'native',              &
     1496                                  MPI_INFO_NULL, ierr )
     1497          CALL MPI_FILE_READ_ALL( fh, array_3di4, SIZE( array_3di4 ), MPI_INTEGER, status, ierr )
     1498       ENDIF
     1499       CALL sm_io%sm_node_barrier()
     1500#else
     1501       CALL posix_lseek( fh, array_position )
     1502       CALL posix_read(fh, array_3di4, SIZE( array_3di4 ) )
     1503#endif
     1504       IF ( include_total_domain_boundaries )  THEN
     1505          DO  i = iog%nxl, iog%nxr
     1506             data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3di4(:,i,iog%nys:iog%nyn)
     1507          ENDDO
     1508       ELSE
     1509          DO  i = nxl, nxr
     1510             data(:,nys:nyn,i) = array_3di4(:,i,nys:nyn)
     1511          ENDDO
     1512       ENDIF
     1513
     1514    ELSE
     1515
     1516       message_string = '3d-INTEGER*4 array "' // TRIM( name ) // '" not found in restart file'
     1517       CALL message( 'rrd_mpi_io_int4_3d', 'PA0722', 3, 2, 0, 6, 0 )
     1518
     1519    ENDIF
     1520
     1521
     1522 END SUBROUTINE rrd_mpi_io_int4_3d
     1523
     1524
     1525
     1526!--------------------------------------------------------------------------------------------------!
     1527! Description:
     1528! ------------
     1529!> Read 3d-INTEGER*8 array with MPI-IO
     1530!--------------------------------------------------------------------------------------------------!
     1531 SUBROUTINE rrd_mpi_io_int8_3d( name, data )
     1532
     1533    IMPLICIT NONE
     1534
     1535    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     1536
     1537    INTEGER(iwp)                       ::  i       !<
     1538
     1539#if defined( __parallel )
     1540    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     1541#endif
     1542
     1543    LOGICAL                            ::  found   !<
     1544
     1545    INTEGER(idp), INTENT(INOUT), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data  !<
     1546
     1547
     1548    found = .FALSE.
     1549    data  = -1.0
     1550
     1551    DO  i = 1, tgh%nr_arrays
     1552       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
     1553          array_position = array_offset(i)
     1554          found = .TRUE.
     1555          EXIT
     1556       ENDIF
     1557    ENDDO
     1558
     1559    IF ( found )  THEN
     1560
     1561#if defined( __parallel )
     1562          CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited # of cores is inactive
     1563          IF( sm_io%iam_io_pe )  THEN
     1564             CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER8, ft_3di8, 'native', MPI_INFO_NULL, &
     1565                                     ierr )
     1566             CALL MPI_FILE_READ_ALL( fh, array_3di8, SIZE( array_3di8 ), MPI_INTEGER8, status, ierr )
     1567          ENDIF
     1568          CALL sm_io%sm_node_barrier()
     1569#else
     1570          CALL posix_lseek( fh, array_position )
     1571          CALL posix_read(fh, array_3di8, SIZE( array_3di8 ) )
     1572#endif
     1573          IF ( include_total_domain_boundaries )  THEN
     1574             DO  i = iog%nxl, iog%nxr
     1575                data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp) = array_3di8(:,i,iog%nys:iog%nyn)
     1576             ENDDO
     1577          ELSE
     1578             DO  i = nxl, nxr
     1579                data(:,nys:nyn,i) = array_3di8(:,i,nys:nyn)
     1580             ENDDO
     1581          ENDIF
     1582
     1583    ELSE
     1584
     1585       message_string = '3d-INTEGER*8 array "' // TRIM( name ) // '" not found in restart file'
     1586       CALL message( 'rrd_mpi_io_int8_3d', 'PA0722', 3, 2, 0, 6, 0 )
     1587
     1588    ENDIF
     1589
     1590
     1591 END SUBROUTINE rrd_mpi_io_int8_3d
     1592
     1593
     1594
     1595!--------------------------------------------------------------------------------------------------!
     1596! Description:
     1597! ------------
    14171598!> Read 2d-REAL array with MPI-IO
    14181599!--------------------------------------------------------------------------------------------------!
    … …  
    19492130! Description:
    19502131! ------------
     2132!> Write 3d-INTEGER*4 array with MPI-IO
     2133!--------------------------------------------------------------------------------------------------!
     2134 SUBROUTINE wrd_mpi_io_int4_3d( name, data )
     2135
     2136    IMPLICIT NONE
     2137
     2138    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     2139
     2140    INTEGER(iwp)                       ::  i       !<
     2141#if defined( __parallel )
     2142    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     2143#endif
     2144    INTEGER(isp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data !<
     2145
     2146
     2147    IF ( header_array_index == max_nr_arrays )  THEN
     2148       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
     2149    ENDIF
     2150
     2151    array_names(header_array_index)  = name
     2152    array_offset(header_array_index) = array_position
     2153    header_array_index = header_array_index + 1
     2154
     2155    IF ( include_total_domain_boundaries )  THEN
     2156!
     2157!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
     2158!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
     2159!--    index order of the array in the same way, i.e. the first dimension should be along x and the
     2160!--    second along y. For this reason, the original PALM data need to be swaped.
     2161       DO  i = iog%nxl, iog%nxr
     2162          array_3di4(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
     2163       ENDDO
     2164
     2165    ELSE
     2166!
     2167!--    Prepare output of 3d-REAL-array without ghost layers
     2168       DO  i = nxl, nxr
     2169           array_3di4(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
     2170       ENDDO
     2171
     2172    ENDIF
     2173#if defined( __parallel )
     2174    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
     2175    IF ( sm_io%iam_io_pe )  THEN
     2176       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER, ft_3di4, 'native', MPI_INFO_NULL, ierr )
     2177       CALL MPI_FILE_WRITE_ALL( fh, array_3di4, SIZE( array_3di4 ), MPI_INTEGER, status, ierr )
     2178    ENDIF
     2179    CALL sm_io%sm_node_barrier()
     2180#else
     2181    CALL posix_lseek( fh, array_position )
     2182    CALL posix_write( fh, array_3di4, SIZE( array_3di4 ) )
     2183#endif
     2184!
     2185!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
     2186!-- Maybe a compiler problem.
     2187    array_position = array_position + INT(     (nz+2), KIND = rd_offset_kind ) *                   &
     2188                                      INT( (iog%ny+1), KIND = rd_offset_kind ) *                   &
     2189                                      INT( (iog%nx+1), KIND = rd_offset_kind ) * isp
     2190
     2191    write(9,*) 'array_position int4_3d ',trim(name),' ',array_position
     2192
     2193 END SUBROUTINE wrd_mpi_io_int4_3d
     2194
     2195
     2196
     2197!--------------------------------------------------------------------------------------------------!
     2198! Description:
     2199! ------------
     2200!> Write 3d-INTEGER*8 array with MPI-IO
     2201!--------------------------------------------------------------------------------------------------!
     2202 SUBROUTINE wrd_mpi_io_int8_3d( name, data )
     2203
     2204    IMPLICIT NONE
     2205
     2206    CHARACTER(LEN=*), INTENT(IN)       ::  name    !<
     2207
     2208    INTEGER(iwp)                       ::  i       !<
     2209#if defined( __parallel )
     2210    INTEGER, DIMENSION(rd_status_size) ::  status  !<
     2211#endif
     2212    INTEGER(idp), INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  data !<
     2213
     2214
     2215    IF ( header_array_index == max_nr_arrays )  THEN
     2216       STOP '+++ maximum number of 2d/3d-array entries in restart file header exceeded'
     2217    ENDIF
     2218
     2219    array_names(header_array_index)  = name
     2220    array_offset(header_array_index) = array_position
     2221    header_array_index = header_array_index + 1
     2222
     2223    IF ( include_total_domain_boundaries )  THEN
     2224!
     2225!--    Prepare output of 3d-REAL-array with ghost layers. In the virtual PE grid, the first
     2226!--    dimension is PEs along x, and the second along y. For MPI-IO it is recommended to have the
     2227!--    index order of the array in the same way, i.e. the first dimension should be along x and the
     2228!--    second along y. For this reason, the original PALM data need to be swaped.
     2229       DO  i = iog%nxl, iog%nxr
     2230          array_3di8(:,i,iog%nys:iog%nyn) = data(:,iog%nys-nbgp:iog%nyn-nbgp,i-nbgp)
     2231       ENDDO
     2232
     2233    ELSE
     2234!
     2235!--    Prepare output of 3d-REAL-array without ghost layers
     2236       DO  i = nxl, nxr
     2237           array_3di8(:,i,iog%nys:iog%nyn) = data(:,nys:nyn,i)
     2238       ENDDO
     2239
     2240    ENDIF
     2241#if defined( __parallel )
     2242    CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
     2243    IF ( sm_io%iam_io_pe )  THEN
     2244       CALL MPI_FILE_SET_VIEW( fh, array_position, MPI_INTEGER8, ft_3di8, 'native', MPI_INFO_NULL, ierr )
     2245       CALL MPI_FILE_WRITE_ALL( fh, array_3di8, SIZE( array_3di8 ), MPI_INTEGER8, status, ierr )
     2246    ENDIF
     2247    CALL sm_io%sm_node_barrier()
     2248#else
     2249    CALL posix_lseek( fh, array_position )
     2250    CALL posix_write( fh, array_3di8, SIZE( array_3di8 ) )
     2251#endif
     2252!
     2253!-- Type conversion required, otherwise right hand side brackets are calculated assuming 4 byte INT.
     2254!-- Maybe a compiler problem.
     2255    array_position = array_position + INT(     (nz+2), KIND = rd_offset_kind ) *                   &
     2256                                      INT( (iog%ny+1), KIND = rd_offset_kind ) *                   &
     2257                                      INT( (iog%nx+1), KIND = rd_offset_kind ) * dp
     2258
     2259    write(9,*) 'array_position int8_3d ',trim(name),' ',array_position
     2260
     2261 END SUBROUTINE wrd_mpi_io_int8_3d
     2262
     2263
     2264
     2265!--------------------------------------------------------------------------------------------------!
     2266! Description:
     2267! ------------
    19512268!> Write 3d-REAL array with MPI-IO
    19522269!--------------------------------------------------------------------------------------------------!
    … …  
    20072324                                      INT( (iog%ny+1), KIND = rd_offset_kind ) *                   &
    20082325                                      INT( (iog%nx+1), KIND = rd_offset_kind ) * wp
     2326
     2327    write(9,*) 'array_position real3d ',trim(name),' ',array_position
    20092328
    20102329 END SUBROUTINE wrd_mpi_io_real_3d
    … …  
    25762895! Description:
    25772896! ------------
     2897!> Read particle data array with MPI-IO.
     2898!--------------------------------------------------------------------------------------------------!
     2899 SUBROUTINE rrd_mpi_io_particles( name, prt_global_index )
     2900
     2901    IMPLICIT NONE
     2902
     2903    CHARACTER(LEN=*), INTENT(IN)                       ::  name            !<
     2904    INTEGER(idp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  prt_global_index      !<
     2905
     2906    INTEGER(iwp)                       ::  array_size      !<
     2907    INTEGER(iwp)                       ::  byte_column     !<
     2908    INTEGER(iwp)                       ::  i               !<
     2909    INTEGER(iwp)                       ::  ind             !<
     2910    INTEGER(iwp)                       ::  j               !<
     2911    INTEGER(iwp)                       ::  k               !<
     2912    INTEGER(iwp)                       ::  n               !<
     2913    INTEGER(iwp)                       ::  particle_size   !<
     2914
     2915    INTEGER(KIND=rd_offset_kind)       ::  disp            !<
     2916    INTEGER(KIND=rd_offset_kind)       ::  offset          !<
     2917    INTEGER(KIND=rd_offset_kind)       ::  prt_nr_bytes    !<
     2918
     2919    LOGICAL                            ::  found           !<
     2920
     2921    REAL(dp)                           :: rr               !< there is no data type INTEGER*8 in MPI
     2922    REAL(dp)                           :: rs               !< use REAL*8 to compute max offset
     2923
     2924    TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET :: prt_data   !<
     2925
     2926#if defined( __parallel )
     2927    INTEGER, DIMENSION(rd_status_size) ::  status          !<
     2928#else
     2929    TYPE(C_PTR)                        ::  buf
     2930#endif
     2931
     2932    found = .FALSE.
     2933
     2934    DO  i = 1, tgh%nr_arrays
     2935       IF ( TRIM(array_names(i)) == TRIM( name ) )  THEN
     2936          array_position = array_offset(i)
     2937          found = .TRUE.
     2938          EXIT
     2939       ENDIF
     2940    ENDDO
     2941
     2942    IF ( found )  THEN
     2943
     2944       offset = 0
     2945
     2946       particle_size = STORAGE_SIZE(zero_particle) / 8  ! 8 here means number of bits/byte, NOT wp
     2947
     2948       array_size = 0
     2949       DO  i = nxl, nxr
     2950          DO  j = nys, nyn
     2951             array_size = MAX( array_size, SUM(prt_count(:,j,i)) )
     2952          ENDDO
     2953       ENDDO
     2954
     2955       write(9,*) 'particle_size_read ',particle_size,array_size,array_position,sum(prt_global_index)
     2956
     2957       ALLOCATE( prt_data(MAX(array_size,1)) )
     2958
     2959!
     2960!--    Write columns of particle
     2961#if defined( __parallel )
     2962       CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
     2963#endif
     2964       prt_nr_bytes = 0
     2965       DO  i = nxl, nxr
     2966          DO  j = nys, nyn
     2967             disp         = array_position + prt_global_index(nzb,j,i) * particle_size
     2968             byte_column  = SUM( prt_count(:,j,i) ) * particle_size
     2969             prt_nr_bytes = MAX( disp+byte_column, prt_nr_bytes )
     2970
     2971#if defined( __parallel )
     2972             CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
     2973             IF ( byte_column > 0 )  THEN
     2974                CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
     2975                CALL MPI_FILE_READ( fh, prt_data, byte_column, MPI_BYTE, status, ierr )
     2976             ENDIF
     2977             CALL sm_io%sm_node_barrier()
     2978#else
     2979             buf = C_LOC(prt_data)     ! use C_PTR to avaid another overlay in posix interface
     2980             CALL posix_lseek( fh, disp )
     2981             CALL posix_read( fh, buf, byte_column )
     2982#endif
     2983             ind = 1
     2984             DO  k = nzb, nzt+1
     2985                DO  n = 1, prt_count(k,j,i)
     2986                   grid_particles(k,j,i)%particles(n) = prt_data(ind)
     2987                   ind = ind+1
     2988                ENDDO
     2989             ENDDO
     2990
     2991          ENDDO
     2992       ENDDO
     2993
     2994#if defined( __parallel )
     2995       rs = prt_nr_bytes
     2996       CALL MPI_ALLREDUCE( rs, rr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm2d, ierr )
     2997       prt_nr_bytes = rr
     2998#else
     2999       rr = rs
     3000#endif
     3001       array_position = prt_nr_bytes
     3002
     3003       write(9,*) 'array_position after particle read ',array_position,prt_nr_bytes,rs
     3004
     3005       DEALLOCATE( prt_data )
     3006
     3007    ENDIF
     3008
     3009 END SUBROUTINE rrd_mpi_io_particles
     3010
     3011
     3012
     3013!--------------------------------------------------------------------------------------------------!
     3014! Description:
     3015! ------------
    25783016!> Read 1d-REAL surface data array with MPI-IO.
    25793017!--------------------------------------------------------------------------------------------------!
    … …  
    25993037#if defined( __parallel )
    26003038    INTEGER, DIMENSION(rd_status_size)  ::  status   !<
     3039#else
     3040    TYPE(C_PTR)                         ::  buf      !<
    26013041#endif
    26023042
    26033043    LOGICAL                             ::  found    !<
    26043044
    2605     REAL(wp), INTENT(OUT), DIMENSION(:) ::  data     !<
     3045    REAL(wp), INTENT(OUT), DIMENSION(:), TARGET ::  data  !<
    26063046
    26073047
    … …  
    26773117                      ierr )
    26783118#else
     3119!
     3120!--                Use C_PTR here, because posix read does not work with indexed array
     3121                   buf = C_LOC( data(m_start_index(j_f,i_f)) )
    26793122                   CALL posix_lseek( fh, disp_f )
    2680                    CALL posix_read( fh, data(m_start_index(j_f,i_f):), nr_bytes_f )
     3123                   CALL posix_read( fh, buf, nr_bytes_f )
    26813124#endif
    26823125                   disp_f     = disp
    … …  
    26973140
    26983141    ENDIF
    2699 
    2700 !    IF ( lo_first_index == 1 )  THEN
    2701 !       IF ( debug_level >= 2 .AND. nr_val > 0 )  WRITE(9,*)  'r_surf_1 ', TRIM( name ), ' ', nr_val, SUM( data(1:nr_val) )
    2702 !    ELSE
    2703 !       IF ( debug_level >= 2 .AND. nr_val > 0 )  WRITE(9,*)  'r_surf_next ', TRIM( name ), ' ', &
    2704 !                                                             lo_first_index,nr_val, SUM( data(1:nr_val) )
    2705 !    ENDIF
    2706 
    27073142
    27083143 CONTAINS
    … …  
    29283363! Description:
    29293364! ------------
     3365!> Write particle data with MPI-IO.
     3366!--------------------------------------------------------------------------------------------------!
     3367 SUBROUTINE wrd_mpi_io_particles( name, prt_global_index )
     3368
     3369    IMPLICIT NONE
     3370
     3371    CHARACTER(LEN=*), INTENT(IN)                       ::  name            !<
     3372    INTEGER(idp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ::  prt_global_index      !<
     3373
     3374    INTEGER(iwp)                       ::  array_size      !<
     3375    INTEGER(iwp)                       ::  byte_column     !<
     3376    INTEGER(iwp)                       ::  i               !<
     3377    INTEGER(iwp)                       ::  ind             !<
     3378    INTEGER(iwp)                       ::  j               !<
     3379    INTEGER(iwp)                       ::  k               !<
     3380    INTEGER(iwp)                       ::  n               !<
     3381    INTEGER(iwp)                       ::  particle_size   !<
     3382
     3383    INTEGER(KIND=rd_offset_kind)       ::  disp            !<
     3384    INTEGER(KIND=rd_offset_kind)       ::  offset          !<
     3385    INTEGER(KIND=rd_offset_kind)       ::  prt_nr_bytes    !<
     3386
     3387    REAL(dp)                           :: rs               !< use REAL*8 to compute max offset
     3388    REAL(dp)                           :: rr               !< there is no data type INTEGER*8 in MPI
     3389
     3390
     3391    TYPE(particle_type), DIMENSION(:), ALLOCATABLE, TARGET ::  prt_data   !<
     3392
     3393#if defined( __parallel )
     3394    INTEGER, DIMENSION(rd_status_size) ::  status          !<
     3395#else
     3396    TYPE(C_PTR)                        ::  buf
     3397#endif
     3398
     3399
     3400    offset = 0
     3401
     3402    array_names(header_array_index)  = TRIM(name)
     3403    array_offset(header_array_index) = array_position
     3404    header_array_index = header_array_index + 1
     3405
     3406    particle_size = STORAGE_SIZE( zero_particle ) / 8
     3407
     3408    array_size = 0
     3409    DO  i = nxl, nxr
     3410      DO  j = nys, nyn
     3411         array_size = MAX( array_size, SUM(prt_count(:,j,i)) )
     3412       ENDDO
     3413    ENDDO
     3414
     3415    ALLOCATE( prt_data(MAX(array_size,1)) )
     3416
     3417!
     3418!-- Write columns of particles.
     3419!-- The particles of a column are stored sequentially in the first dimension of the particles array.
     3420!-- Store only the particles of one cell would be possible with this setup, but the I/O portions
     3421!-- for a maximum of 100 particles are not big enough.
     3422#if defined( __parallel )
     3423    CALL MPI_FILE_SET_VIEW( fh, offset, MPI_BYTE, MPI_BYTE, 'native', MPI_INFO_NULL, ierr )
     3424#endif
     3425    prt_nr_bytes = 0
     3426    DO  i = nxl, nxr
     3427       DO  j = nys, nyn
     3428          disp         = array_position + prt_global_index(nzb,j,i) * particle_size
     3429          byte_column  = SUM( prt_count(:,j,i) ) * particle_size
     3430          prt_nr_bytes = MAX( disp+byte_column, prt_nr_bytes )
     3431
     3432          ind = 1
     3433          DO  k = nzb, nzt+1
     3434             DO  n = 1, prt_count(k,j,i)
     3435                prt_data(ind) = grid_particles(k,j,i)%particles(n)
     3436                ind = ind+1
     3437             ENDDO
     3438          ENDDO
     3439
     3440#if defined( __parallel )
     3441          CALL sm_io%sm_node_barrier()  ! Has no effect if I/O on limited number of cores is inactive
     3442          IF ( byte_column > 0 )  THEN
     3443             CALL MPI_FILE_SEEK( fh, disp, MPI_SEEK_SET, ierr )
     3444             CALL MPI_FILE_WRITE( fh, prt_data, byte_column, MPI_BYTE, status, ierr )
     3445          ENDIF
     3446          CALL sm_io%sm_node_barrier()
     3447#else
     3448          buf = C_LOC(prt_data)  ! use C_PTR to avoid another overlay in posix interface
     3449          CALL posix_lseek( fh, disp )
     3450          CALL posix_write( fh, buf, byte_column )
     3451#endif
     3452       ENDDO
     3453    ENDDO
     3454
     3455#if defined( __parallel )
     3456    rs = prt_nr_bytes
     3457    CALL MPI_ALLREDUCE( rs, rr, 1, MPI_DOUBLE_PRECISION, MPI_MAX, comm2d, ierr )
     3458    prt_nr_bytes = rr
     3459#else
     3460    rr = rs
     3461#endif
     3462    array_position = prt_nr_bytes
     3463
     3464    write(9,*) 'array_position after particle ',array_position,prt_nr_bytes,rs
     3465
     3466    DEALLOCATE( prt_data )
     3467
     3468 END SUBROUTINE wrd_mpi_io_particles
     3469
     3470
     3471
     3472!--------------------------------------------------------------------------------------------------!
     3473! Description:
     3474! ------------
    29303475!> Write 1d-REAL surface data array with MPI-IO.
    29313476!--------------------------------------------------------------------------------------------------!
    … …  
    30433588
    30443589
     3590
    30453591!--------------------------------------------------------------------------------------------------!
    30463592! Description:
    … …  
    32073753    ENDIF
    32083754
     3755    fh = -1
     3756
    32093757    restart_file_size = array_position / ( 1024.0_dp * 1024.0_dp )
    32103758
    32113759 END SUBROUTINE rd_mpi_io_close
     3760
     3761
     3762
     3763 FUNCTION rd_mpi_io_check_open()  RESULT ( isopen )
     3764
     3765    LOGICAL ::  isopen
     3766
     3767    isopen = ( fh /= -1 )
     3768
     3769 END FUNCTION rd_mpi_io_check_open
    32123770
    32133771
    … …  
    35674125! Description:
    35684126! ------------
     4127!> This subroutine creates file types to access 3d-INTEGER*4 arrays and 3d-INTEGER*8 arrays
     4128!> distributed in blocks among processes to a single file that contains the global arrays.
     4129!> These types are only used for particle data.
     4130!--------------------------------------------------------------------------------------------------!
     4131  SUBROUTINE rd_mpi_io_particle_filetypes
     4132
     4133    IMPLICIT NONE
     4134
     4135    INTEGER, DIMENSION(3) ::  dims3   !<
     4136    INTEGER, DIMENSION(3) ::  lize3   !<
     4137    INTEGER, DIMENSION(3) ::  start3  !<
     4138
     4139    TYPE(domain_decomposition_grid_features) ::  save_io_grid  !< temporary variable to store grid settings
     4140
     4141!
     4142!-- MPI_INTEGER8 is not standard MPI, but is supported on most MPI distibutions
     4143!-- If not suppported, a workaround could be enabled with the following preprocessor directive
     4144!#if defined( __NO_INTEGER8)
     4145!    CALL MPI_TYPE_CONTIGUOUS( 2, MPI_INTEGER, MPI_INTEGER8, ierr )
     4146!    CALL MPI_TYPE_COMMIT( MPI_INTEGER8, ierr )
     4147!#endif
     4148
     4149    IF ( sm_io%is_sm_active() )  THEN
     4150       save_io_grid = sm_io%io_grid
     4151    ENDIF
     4152
     4153    IF( .NOT. pe_active_for_read )  RETURN
     4154
     4155!
     4156!-- Decision, if storage with or without ghost layers.
     4157!-- Please note that the indexing of the global array always starts at 0, even in Fortran.
     4158!-- Therefore the PE local indices have to be shifted by nbgp in the case with ghost layers.
     4159    IF ( include_total_domain_boundaries )  THEN
     4160
     4161       iog%nxl = nxl + nbgp
     4162       iog%nxr = nxr + nbgp
     4163       iog%nys = nys + nbgp
     4164       iog%nyn = nyn + nbgp
     4165       iog%nnx = nnx
     4166       iog%nny = nny
     4167       iog%nx  = nx + 2 * nbgp
     4168       iog%ny  = ny + 2 * nbgp
     4169       IF ( myidx == 0 )  THEN
     4170          iog%nxl = iog%nxl - nbgp
     4171          iog%nnx = iog%nnx + nbgp
     4172       ENDIF
     4173       IF ( myidx == npex-1  .OR.  npex == -1 )  THEN   ! npex == 1 if -D__parallel not set
     4174          iog%nxr = iog%nxr + nbgp
     4175          iog%nnx = iog%nnx + nbgp
     4176       ENDIF
     4177       IF ( myidy == 0 )  THEN
     4178          iog%nys = iog%nys - nbgp
     4179          iog%nny = iog%nny + nbgp
     4180       ENDIF
     4181       IF ( myidy == npey-1  .OR.  npey == -1 )  THEN   ! npey == 1 if -D__parallel not set
     4182          iog%nyn = iog%nyn + nbgp
     4183          iog%nny = iog%nny + nbgp
     4184       ENDIF
     4185
     4186       CALL sm_io%sm_adjust_outer_boundary()
     4187
     4188    ELSE
     4189
     4190       iog%nxl = nxl
     4191       iog%nxr = nxr
     4192       iog%nys = nys
     4193       iog%nyn = nyn
     4194       iog%nnx = nnx
     4195       iog%nny = nny
     4196       iog%nx  = nx
     4197       iog%ny  = ny
     4198
     4199    ENDIF
     4200
     4201    IF ( sm_io%is_sm_active() )  THEN
     4202#if defined( __parallel )
     4203       CALL sm_io%sm_allocate_shared( array_3di4, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,&
     4204                                      sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3di4 )
     4205       CALL sm_io%sm_allocate_shared( array_3di8, nzb, nzt+1, sm_io%io_grid%nxl, sm_io%io_grid%nxr,&
     4206                                      sm_io%io_grid%nys, sm_io%io_grid%nyn, win_3di8 )
     4207#endif
     4208    ELSE
     4209       ALLOCATE( array_3di4(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
     4210       ALLOCATE( array_3di8(nzb:nzt+1,iog%nxl:iog%nxr,iog%nys:iog%nyn) )
     4211
     4212       sm_io%io_grid = iog
     4213    ENDIF
     4214
     4215!
     4216!-- Create filetype for 3d INTEGER array
     4217    dims3(1)  = nz + 2
     4218    dims3(2)  = iog%nx + 1
     4219    dims3(3)  = iog%ny + 1
     4220
     4221    lize3(1)  = dims3(1)
     4222    lize3(2)  = sm_io%io_grid%nnx
     4223    lize3(3)  = sm_io%io_grid%nny
     4224
     4225    start3(1) = nzb
     4226    start3(2) = sm_io%io_grid%nxl
     4227    start3(3) = sm_io%io_grid%nys
     4228
     4229#if defined( __parallel )
     4230    IF ( sm_io%iam_io_pe )  THEN
     4231       CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_INTEGER,     &
     4232                                      ft_3di4, ierr )
     4233       CALL MPI_TYPE_COMMIT( ft_3di4, ierr )
     4234
     4235       CALL MPI_TYPE_CREATE_SUBARRAY( 3, dims3, lize3, start3, MPI_ORDER_FORTRAN, MPI_INTEGER8,    &
     4236                                      ft_3di8, ierr )
     4237       CALL MPI_TYPE_COMMIT( ft_3di8, ierr )
     4238    ENDIF
     4239#endif
     4240
     4241 END SUBROUTINE rd_mpi_io_particle_filetypes
     4242
     4243
     4244
     4245!--------------------------------------------------------------------------------------------------!
     4246! Description:
     4247! ------------
    35694248!> This subroutine creates file types to access 3d-soil arrays distributed in blocks among processes
    35704249!> to a single file that contains the global arrays. It is not required for the serial mode.
    … …  
    36554334    ENDIF
    36564335
     4336!
     4337!-- Free last particle filetypes
     4338    IF ( sm_io%iam_io_pe .AND. ft_3di4 /= -1 )  THEN
     4339       CALL MPI_TYPE_FREE( ft_3di4, ierr )
     4340       CALL MPI_TYPE_FREE( ft_3di8, ierr )
     4341    ENDIF
     4342
     4343    IF ( sm_io%is_sm_active() .AND.  win_3di4 /= -1 )  THEN
     4344       CALL sm_io%sm_free_shared( win_3di4 )
     4345       CALL sm_io%sm_free_shared( win_3di8 )
     4346    ENDIF
     4347
    36574348    ft_surf  = -1
    36584349    win_surf = -1
    36594350#else
    3660     DEALLOCATE( array_2d, array_2di, array_3d )
     4351    IF ( ASSOCIATED(array_2d)   )  DEALLOCATE( array_2d )
     4352    IF ( ASSOCIATED(array_2di)  )  DEALLOCATE( array_2di )
     4353    IF ( ASSOCIATED(array_3d)   )  DEALLOCATE( array_3d )
     4354    IF ( ASSOCIATED(array_3di4) )  DEALLOCATE( array_3di4 )
     4355    IF ( ASSOCIATED(array_3di8) )  DEALLOCATE( array_3di8 )
    36614356#endif
    36624357
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |