Ignore:
Timestamp:
Nov 26, 2019 6:49:59 PM (17 months ago)
Author:
suehring
Message:

Synthetic turbulence generator: Computation of velocity seeds optimized. This implies that random numbers are computed now using the parallel random number generator. Random number are now only computed and normalized locally, while distributed over all mpi ranks afterwards, instead of computing random numbers on a global array. urther, the number of calls for the time-consuming velocity-seed generation is reduced - now the left and right, as well as the north and south boundary share the same velocity-seed matrices.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

    r4182 r4309  
    2525! -----------------
    2626! $Id$
     27! Computation of velocity seeds optimized. This implies that random numbers
     28! are computed now using the parallel random number generator. Random numbers
     29! are now only computed and normalized locally, while distributed over all 
     30! mpi ranks afterwards, instead of computing random numbers on a global array.
     31! Further, the number of calls for the time-consuming velocity-seed generation
     32! is reduced - now the left and right, as well as the north and south boundary
     33! share the same velocity-seed matrices.
     34!
     35! 4182 2019-08-22 15:20:23Z scharf
    2736! Corrected "Former revisions" section
    2837!
     
    154163               neutral,                                                        &
    155164               num_mean_inflow_profiles,                                       &
     165               random_generator,                                               &
    156166               rans_mode,                                                      &
    157167               restart_string,                                                 &
     
    206216    USE pmc_interface,                                                         &
    207217        ONLY : rans_mode_parent
     218       
     219   USE random_generator_parallel,                                              &
     220        ONLY:  random_dummy,                                                   &
     221               random_number_parallel,                                         &
     222               random_seed_parallel,                                           &
     223               seq_random_array
    208224
    209225    USE transpose_indices,                                                     &
     
    227243    INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
    228244    INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
    229     INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
    230     INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
    231     INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
    232     INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
    233245    INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
    234246    INTEGER(iwp) ::  mergp              !< merg + nbgp
     
    237249    INTEGER(iwp) ::  nzb_y_stg          !< lower bound of z coordinate (required for transposing z on PEs along y)
    238250    INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
     251    INTEGER(iwp) ::  stg_type_xz        !< MPI type for full z range
     252    INTEGER(iwp) ::  stg_type_xz_small  !< MPI type for small z range
     253    INTEGER(iwp) ::  stg_type_yz        !< MPI type for full z range
     254    INTEGER(iwp) ::  stg_type_yz_small  !< MPI type for small z range
    239255
    240256    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections,
     
    256272    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwy            !< length scale of w in y direction (in gp)
    257273    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwz            !< length scale of w in z direction (in gp)
    258     INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
    259274
    260275    REAL(wp) ::  blend                    !< value to create gradually and smooth blending of Reynolds stress and length
     
    476491          CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
    477492       ENDIF
     493!
     494!--    Synthetic turbulence generator requires the parallel random generator
     495       IF ( random_generator /= 'random-parallel' )  THEN
     496          message_string = 'Using synthetic turbulence generator ' //          &
     497                           'requires random_generator = random-parallel.'
     498          CALL message( 'stg_check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
     499       ENDIF
    478500
    479501    ENDIF
     
    540562    INTEGER(iwp) :: newtype                  !< dummy MPI type
    541563    INTEGER(iwp) :: realsize                 !< size of REAL variables
    542     INTEGER(iwp) :: nseed                    !< dimension of random number seed
    543     INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
    544    
     564
    545565    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz_l = 0 !< number of non-topography grid points at xz-cross-section on subdomain
    546566    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_yz_l = 0 !< number of non-topography grid points at yz-cross-section on subdomain
     
    648668
    649669#endif
    650 !
    651 !-- Define seed of random number
    652     CALL RANDOM_SEED()
    653     CALL RANDOM_SEED( size=nseed )
    654     ALLOCATE( seed(1:nseed) )
    655     DO  j = 1, nseed
    656        seed(j) = startseed + j
    657     ENDDO
    658     CALL RANDOM_SEED(put=seed)
    659670!
    660671!-- Allocate required arrays.
     
    12091220!-- Initial value of fu, fv, fw
    12101221    IF ( time_since_reference_point == 0.0_wp .AND. .NOT. velocity_seed_initialized )  THEN
    1211        CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
    1212        CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
    1213        CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
    1214 
    1215        IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
    1216                                 .AND.  .NOT.  rans_mode ) )  THEN
    1217 !
    1218 !--       Generate turbulence at right boundary
    1219           CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_right )
    1220           CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_right )
    1221           CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_right )
    1222 !
    1223 !--       Generate turbulence at north boundary
    1224           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_north )
    1225           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_north )
    1226           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_north )
    1227 !
    1228 !--       Generate turbulence at south boundary
    1229           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz, id_stg_south )
    1230           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz, id_stg_south )
    1231           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz, id_stg_south )
     1222!
     1223!--    Generate turbulence at the left and right boundary. Random numbers
     1224!--    for the yz-planes at the left/right boundary are generated by the
     1225!--    left-sided mpi ranks only. After random numbers are calculated, they
     1226!--    are distributed to all other mpi ranks in the model, so that the
     1227!--    velocity seed matrices are available on all mpi ranks (i.e. also on the
     1228!--    right-sided boundary mpi ranks). In case of offline nesting, this implies,
     1229!--    that the left- and the right-sided lateral boundary have the same initial
     1230!--    seeds.
     1231!--    Note, in case of inflow from the right only, only turbulence at the left
     1232!--    boundary is required.
     1233       IF ( .NOT. ( nesting_offline  .OR.                                      &
     1234                  ( child_domain .AND.  rans_mode_parent                       &
     1235                                 .AND.  .NOT.  rans_mode ) ) )  THEN
     1236
     1237          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
     1238          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
     1239          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz, id_stg_left )
     1240       ELSE
     1241          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz,               &
     1242                                     id_stg_left, id_stg_right )
     1243          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz,               &
     1244                                     id_stg_left, id_stg_right )
     1245          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw_yz,               &
     1246                                     id_stg_left, id_stg_right )
     1247!
     1248!--       Generate turbulence at the south and north boundary. Random numbers
     1249!--       for the xz-planes at the south/north boundary are generated by the
     1250!--       south-sided mpi ranks only. Please see also comment above.
     1251          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fu_xz,               &
     1252                                     id_stg_south, id_stg_north )
     1253          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fv_xz,               &
     1254                                     id_stg_south, id_stg_north )
     1255          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fw_xz,               &
     1256                                     id_stg_south, id_stg_north )
    12321257       ENDIF
    12331258       velocity_seed_initialized = .TRUE.
    12341259    ENDIF
    1235    
    1236 !
    1237 !-- New set of fu, fv, fw
     1260!
     1261!-- New set of fu, fv, fw. Note, for inflow from the left side only, velocity
     1262!-- seeds are only required at the left boundary, while in case of offline
     1263!-- nesting or RANS-LES nesting velocity seeds are required also at the
     1264!-- right, south and north boundaries.
    12381265    IF ( stg_call )  THEN
    1239 !     
    1240 !--    Generate turbulence at left boundary (required in all applications,
    1241 !--    i.e. offline nesting and turbulent inflow from the left)
    1242        CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
    1243        CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
    1244        CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
    1245        
    1246        IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
    1247                              .AND.  .NOT.  rans_mode ) )  THEN
    1248 !     
    1249 !--       Generate turbulence at right boundary
    1250           CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
    1251           CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
    1252           CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
    1253 !     
    1254 !--       Generate turbulence at north boundary
    1255           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
    1256           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
    1257           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
    1258 !     
    1259 !--       Generate turbulence at south boundary
    1260           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
    1261           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
    1262           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     1266       IF ( .NOT. ( nesting_offline  .OR.                                      &
     1267                  ( child_domain .AND.  rans_mode_parent                       &
     1268                                 .AND.  .NOT.  rans_mode ) ) )  THEN
     1269          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
     1270          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
     1271          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
     1272
     1273       ELSE
     1274          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz,               &
     1275                                     id_stg_left, id_stg_right )
     1276          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz,               &
     1277                                     id_stg_left, id_stg_right )
     1278          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz,               &
     1279                                     id_stg_left, id_stg_right )
     1280
     1281          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz,               &
     1282                                     id_stg_south, id_stg_north )
     1283          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz,               &
     1284                                     id_stg_south, id_stg_north )
     1285          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz,               &
     1286                                     id_stg_south, id_stg_north )
    12631287       ENDIF
    12641288    ENDIF
    12651289   
    12661290!
    1267 !-- Turbulence generation at left and or right boundary
     1291!-- Turbulence generation at left and/or right boundary
    12681292    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
    12691293!
     
    13131337                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) *          &
    13141338                                    MERGE( 1.0_wp, 0.0_wp,                     &
    1315                                            BTEST( wall_flags_0(k,j,i), 1 ) )   
     1339                                           BTEST( wall_flags_0(k,j,i), 1 ) )
    13161340             ENDDO
    13171341          ENDDO
    1318          
     1342
    13191343          IF ( myidx == id_stg_left  )  i = nxl-1
    13201344          IF ( myidx == id_stg_right )  i = nxr+1
     
    13411365          ENDDO
    13421366       ENDIF
    1343 
    13441367!
    13451368!--    Mass flux correction following Kim et al. (2013)
     
    15141537                                         + a22(k) * fv_xz(k,i) ), 3.0_wp ) *   &
    15151538                                    MERGE( 1.0_wp, 0.0_wp,                     &
    1516                                            BTEST( wall_flags_0(k,j,i), 2 ) )   
     1539                                           BTEST( wall_flags_0(k,j,i), 2 ) )
    15171540             ENDDO
    15181541          ENDDO
     
    16481671!> parts are collected to form the full array.
    16491672!------------------------------------------------------------------------------!
    1650  SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
     1673 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id_left, id_right )
     1674 
     1675 USE pegrid
    16511676
    16521677    IMPLICIT NONE
    16531678
    1654     INTEGER(iwp) :: id          !< core ids at respective boundaries
    1655     INTEGER(iwp) :: j           !< loop index in y-direction
    1656     INTEGER(iwp) :: jj          !< loop index in y-direction
    1657     INTEGER(iwp) :: k           !< loop index in z-direction
    1658     INTEGER(iwp) :: kk          !< loop index in z-direction
    1659     INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
     1679    INTEGER(iwp)           :: i           !< grid index x-direction
     1680    INTEGER(iwp)           :: id_left     !< core ids at respective boundaries
     1681    INTEGER(iwp), OPTIONAL :: id_right    !< core ids at respective boundaries
     1682    INTEGER(iwp)           :: j           !< loop index in y-direction
     1683    INTEGER(iwp)           :: jj          !< loop index in y-direction
     1684    INTEGER(iwp)           :: k           !< loop index in z-direction
     1685    INTEGER(iwp)           :: kk          !< loop index in z-direction
     1686    INTEGER(iwp)           :: send_count  !< send count for MPI_GATHERV
    16601687
    16611688    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
     
    16701697    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
    16711698    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
     1699   
    16721700    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
    1673     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
    1674     REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
    1675 
    1676 
    1677 !
    1678 !-- Generate random numbers using a seed generated in stg_init.
     1701    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)             :: f_n     !<  velocity seed
     1702   
     1703    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
     1704    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it_l !< local array of random numbers
     1705
     1706!
     1707!-- Generate random numbers using the parallel random generator.
    16791708!-- The set of random numbers are modified to have an average of 0 and
    1680 !-- unit variance.
     1709!-- unit variance. Note, at the end the random number array must be defined
     1710!-- globally in order to compute the correlation matrices. However,
     1711!-- the calculation and normalization of random numbers is done locally,
     1712!-- while the result is later distributed to all processes. Further,
     1713!-- please note, a set of random numbers is only calculated for the
     1714!-- left boundary, while the right boundary uses the same random numbers
     1715!-- and thus also computes the same correlation matrix.
    16811716    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
     1717    rand_it   = 0.0_wp
    16821718
    16831719    rand_av        = 0.0_wp
    16841720    rand_sigma_inv = 0.0_wp
    16851721    nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
    1686 
    1687     DO  j = 0, ny
    1688        DO  k = nzb, nzt+1
    1689           CALL RANDOM_NUMBER( rand_it(k,j) )
    1690           rand_av = rand_av + rand_it(k,j)
     1722!
     1723!-- Compute and normalize random numbers only on left-boundary ranks.
     1724    IF ( myidx == id_stg_left )  THEN
     1725!
     1726!--    Allocate array for local set of random numbers
     1727       ALLOCATE( rand_it_l(nzb:nzt+1,nys:nyn) )
     1728       rand_it_l = 0.0_wp
     1729
     1730       i = nxl
     1731       DO  j = nys, nyn
     1732!
     1733!--       Put the random seeds at grid point j,i
     1734          CALL random_seed_parallel( put=seq_random_array(:, j, i) )
     1735          DO  k = nzb, nzt+1
     1736             CALL random_number_parallel( random_dummy )
     1737             rand_it_l(k,j) = random_dummy
     1738             rand_av        = rand_av + rand_it_l(k,j)
     1739          ENDDO
     1740!
     1741!--       Get the new random seeds from last call at grid point j,i
     1742          CALL random_seed_parallel( get=seq_random_array(:, j, i) )
    16911743       ENDDO
    1692     ENDDO
    1693 
    1694     rand_av = rand_av * nyz_inv
    1695 
    1696     DO  j = 0, ny
    1697        DO  k = nzb, nzt+1
    1698           rand_it(k,j)   = rand_it(k,j) - rand_av
    1699           rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
     1744#if defined( __parallel )
     1745!
     1746!--    Sum-up the local averages of the random numbers
     1747       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                 &
     1748                           MPI_SUM, comm1dy, ierr )
     1749#endif
     1750       rand_av = rand_av * nyz_inv
     1751!
     1752!--    Obtain zero mean
     1753       rand_it_l = rand_it_l - rand_av
     1754!
     1755!--    Now, compute the variance
     1756       DO  j = nys, nyn
     1757          DO  k = nzb, nzt+1
     1758             rand_sigma_inv = rand_sigma_inv + rand_it_l(k,j)**2
     1759          ENDDO
    17001760       ENDDO
    1701     ENDDO
    1702 
    1703     rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
    1704 
    1705     DO  j = 0, ny
    1706        DO  k = nzb, nzt+1
    1707           rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
    1708        ENDDO
    1709     ENDDO
    1710 
    1711 !
    1712 !-- Periodic fill of random number in space
     1761
     1762#if defined( __parallel )
     1763!
     1764!--    Sum-up the local quadratic averages of the random numbers
     1765       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
     1766                           MPI_SUM, comm1dy, ierr )
     1767#endif
     1768!
     1769!--    Compute standard deviation
     1770       IF ( rand_sigma_inv /= 0.0_wp )  THEN
     1771          rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nyz_inv )
     1772       ELSE
     1773          rand_sigma_inv = 1.0_wp
     1774       ENDIF
     1775!
     1776!--    Normalize with standard deviation to obtain unit variance
     1777       rand_it_l = rand_it_l * rand_sigma_inv
     1778!
     1779!--    Copy local random numbers on the global array
     1780       rand_it(nzb:nzt+1,nys:nyn) = rand_it_l(nzb:nzt+1,nys:nyn)
     1781!
     1782!--    Deallocate local array
     1783       DEALLOCATE( rand_it_l )
     1784!
     1785!--    Now, distribute the final set of random numbers to all mpi ranks located
     1786!--    on the left boundary. Here, an allreduce with sum reduction is sufficient,
     1787!--    or, in the non-parallel case, nothing need to be done at all.
     1788#if defined( __parallel )
     1789       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,   &
     1790                           MPI_SUM, comm1dy, ierr )
     1791#endif
     1792    ENDIF
     1793!
     1794!-- Finally, distribute the set of random numbers (defined on the leftmost-
     1795!-- located mpi ranks) to all other mpi ranks. Here, a allreduce with sum
     1796!-- option is sufficient, because rand_it is zero on all other mpi_ranks.
     1797!-- Note, the reduce operation is only performed with communicator comm1dx,
     1798!-- where only 1 rank within the communicator has non-zero random numbers.
     1799#if defined( __parallel )
     1800    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,      &
     1801                        MPI_SUM, comm1dx, ierr )
     1802#endif
     1803!
     1804!-- Periodic fill of random numbers in space
    17131805    DO  j = 0, ny
    17141806       DO  k = 1, mergp
     
    17481840
    17491841#if defined( __parallel )
    1750     CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,    &
    1751                       f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,  &
    1752                       id, comm1dx, ierr )
     1842!
     1843!-- Gather the velocity seed matrix on left boundary mpi ranks.
     1844    CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,     &
     1845                      f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,   &
     1846                      id_left, comm1dx, ierr )
     1847!
     1848!-- If required, gather the same velocity seed matrix on right boundary
     1849!-- mpi ranks (in offline nesting for example).
     1850    IF ( PRESENT( id_right ) )  THEN
     1851       CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,  &
     1852                         f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,&
     1853                         id_right, comm1dx, ierr )
     1854    ENDIF
    17531855#else
    17541856    f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
     
    17571859
    17581860 END SUBROUTINE stg_generate_seed_yz
    1759 
    1760 
    17611861
    17621862
     
    17701870!> parts are collected to form the full array.
    17711871!------------------------------------------------------------------------------!
    1772  SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
     1872 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north )
    17731873
    17741874    IMPLICIT NONE
    17751875
    1776     INTEGER(iwp) :: id          !< core ids at respective boundaries
    17771876    INTEGER(iwp) :: i           !< loop index in x-direction
     1877    INTEGER(iwp) :: id_north    !< core ids at respective boundaries
     1878    INTEGER(iwp) :: id_south    !< core ids at respective boundaries
    17781879    INTEGER(iwp) :: ii          !< loop index in x-direction
     1880    INTEGER(iwp) :: j           !< grid index y-direction
    17791881    INTEGER(iwp) :: k           !< loop index in z-direction
    17801882    INTEGER(iwp) :: kk          !< loop index in z-direction
     
    17921894    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
    17931895    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
     1896   
    17941897    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
    1795     REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
    1796     REAL(wp), DIMENSION(:,:), ALLOCATABLE        :: rand_it !<  random number
    1797 
    1798 
    1799 !
    1800 !-- Generate random numbers using a seed generated in stg_init.
     1898    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)             :: f_n     !<  velocity seed
     1899
     1900    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
     1901    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it_l !< local array of random numbers
     1902
     1903!
     1904!-- Generate random numbers using the parallel random generator.
    18011905!-- The set of random numbers are modified to have an average of 0 and
    1802 !-- unit variance.
     1906!-- unit variance. Note, at the end the random number array must be defined
     1907!-- globally in order to compute the correlation matrices. However,
     1908!-- the calculation and normalization of random numbers is done locally,
     1909!-- while the result is later distributed to all processes. Further,
     1910!-- please note, a set of random numbers is only calculated for the
     1911!-- left boundary, while the right boundary uses the same random numbers
     1912!-- and thus also computes the same correlation matrix.
    18031913    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
     1914    rand_it   = 0.0_wp
    18041915
    18051916    rand_av        = 0.0_wp
    18061917    rand_sigma_inv = 0.0_wp
    18071918    nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
    1808 
    1809     DO  i = 0, nx
    1810        DO  k = nzb, nzt+1
    1811           CALL RANDOM_NUMBER( rand_it(k,i) )
    1812           rand_av = rand_av + rand_it(k,i)
     1919!
     1920!-- Compute and normalize random numbers only on south-boundary ranks.
     1921    IF ( myidy == id_stg_south )  THEN
     1922!
     1923!--    Allocate array for local set of random numbers
     1924       ALLOCATE( rand_it_l(nzb:nzt+1,nxl:nxr) )
     1925       rand_it_l = 0.0_wp
     1926
     1927       j = nys
     1928       DO  i = nxl, nxr
     1929!
     1930!--       Put the random seeds at grid point j,i
     1931          CALL random_seed_parallel( put=seq_random_array(:, j, i) )
     1932          DO  k = nzb, nzt+1
     1933             CALL random_number_parallel( random_dummy )
     1934             rand_it_l(k,i) = random_dummy
     1935             rand_av = rand_av + rand_it_l(k,i)
     1936          ENDDO
     1937!
     1938!--       Get the new random seeds from last call at grid point j,i
     1939          CALL random_seed_parallel( get=seq_random_array(:, j, i) )
    18131940       ENDDO
    1814     ENDDO
    1815 
    1816     rand_av = rand_av * nxz_inv
    1817 
    1818     DO  i = 0, nx
    1819        DO  k = nzb, nzt+1
    1820           rand_it(k,i)   = rand_it(k,i) - rand_av
    1821           rand_sigma_inv = rand_sigma_inv + rand_it(k,i) ** 2
     1941#if defined( __parallel )
     1942       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                 &
     1943                           MPI_SUM, comm1dx, ierr )
     1944#endif
     1945       rand_av = rand_av * nxz_inv
     1946!
     1947!--    Obtain zero mean
     1948       rand_it_l = rand_it_l - rand_av
     1949!
     1950!--    Now, compute the variance
     1951       DO  i = nxl, nxr
     1952          DO  k = nzb, nzt+1
     1953             rand_sigma_inv = rand_sigma_inv + rand_it_l(k,i)**2
     1954          ENDDO
    18221955       ENDDO
    1823     ENDDO
    1824 
    1825     rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nxz_inv)
    1826 
    1827     DO  i = 0, nx
    1828        DO  k = nzb, nzt+1
    1829           rand_it(k,i) = rand_it(k,i) * rand_sigma_inv
    1830        ENDDO
    1831     ENDDO
    1832 
     1956
     1957#if defined( __parallel )
     1958       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
     1959                           MPI_SUM, comm1dx, ierr )
     1960#endif
     1961!
     1962!--    Compute standard deviation
     1963       IF ( rand_sigma_inv /= 0.0_wp )  THEN
     1964          rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nxz_inv )
     1965       ELSE
     1966          rand_sigma_inv = 1.0_wp
     1967       ENDIF
     1968!
     1969!--    Normalize with standard deviation to obtain unit variance
     1970       rand_it_l = rand_it_l * rand_sigma_inv
     1971!
     1972!--    Copy local random numbers on the global array
     1973       rand_it(nzb:nzt+1,nxl:nxr) = rand_it_l(nzb:nzt+1,nxl:nxr)
     1974!
     1975!--    Deallocate local array
     1976       DEALLOCATE( rand_it_l )
     1977!
     1978!--    Now, distribute the final set of random numbers to all mpi ranks located
     1979!--    on the south boundary. Here, an allreduce with sum reduction is sufficient.
     1980#if defined( __parallel )
     1981       CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,   &
     1982                           MPI_SUM, comm1dx, ierr )
     1983#endif
     1984    ENDIF
     1985!
     1986!-- Finally, distribute the set of random numbers (defined on the 
     1987!-- southmost-located mpi ranks) to all other mpi ranks. Here, a allreduce
     1988!-- with sum option is sufficient, because rand_it is zero on all other
     1989!-- mpi_ranks. Note, the reduce operation is only performed with communicator
     1990!-- comm1dy, where only 1 rank within the communicator has non zero random numbers.
     1991#if defined( __parallel )
     1992    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,      &
     1993                        MPI_SUM, comm1dy, ierr )
     1994#endif
    18331995!
    18341996!-- Periodic fill of random number in space
     
    18452007       ENDDO
    18462008    ENDDO
    1847 
    18482009!
    18492010!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
     
    18702031    IF ( nzt_y_stg == nzt )  send_count = send_count + 1
    18712032
    1872 
    18732033#if defined( __parallel )
    1874     CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,    &
    1875                       f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz,  &
    1876                       id, comm1dy, ierr )
     2034!
     2035!-- Gather the processed velocity seed on south boundary mpi ranks.
     2036    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,   &
     2037                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
     2038                      id_south, comm1dy, ierr )
     2039!
     2040!-- Gather the processed velocity seed on north boundary mpi ranks.
     2041    CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,   &
     2042                      f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
     2043                      id_north, comm1dy, ierr )
    18772044#else
    18782045    f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
    18792046#endif
    1880 
    18812047
    18822048 END SUBROUTINE stg_generate_seed_xz
Note: See TracChangeset for help on using the changeset viewer.