Home


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

Ignore:
Timestamp:
Mar 3, 2020 8:49:28 PM (5 years ago)
Author:
suehring
Message:

Synthetic turbulence: performance optimizations - random numbers only defined and computed locally, option to compute velocity seeds locally without need of global communication; paralell random number generator: new routine to initialize 1D random number arrays; virtual measurements: CPU-log points added

File:
1 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

    r4346 r4438  
    2525! -----------------
    2626! $Id$
     27! Performance optimizations in velocity-seed calculation:
     28!  - random number array is only defined and computed locally (except for
     29!    normalization to zero mean and unit variance)
     30!  - parallel random number generator is applied independent on the 2D random
     31!    numbers in other routines
     32!  - option to decide wheter velocity seeds are computed locally without any
     33!    further communication or are computed by all processes along the
     34!    communicator
     35!
     36! 4346 2019-12-18 11:55:56Z motisi
    2737! Introduction of wall_flags_total_0, which currently sets bits based on static
    2838! topography information used in wall_flags_static_0
    … …  
    184194               time_since_reference_point,                                     &
    185195               turbulent_inflow
    186                
     196
     197    USE cpulog,                                                                &
     198        ONLY:  cpu_log,                                                        &
     199               log_point_s
     200
    187201    USE grid_variables,                                                        &
    188202        ONLY:  ddx,                                                            &
    … …  
    199213               nxl,                                                            &
    200214               nxlu,                                                           &
    201                nxlg,                                                           &
    202215               nxr,                                                            &
    203                nxrg,                                                           &
    204216               ny,                                                             &
    205217               nys,                                                            &
    206218               nysv,                                                           &
    207219               nyn,                                                            &
    208                nyng,                                                           &
    209                nysg,                                                           &
    210220               wall_flags_total_0
    211221
    … …  
    228238               myidy,                                                          &
    229239               pdims
    230        
     240
    231241    USE pmc_interface,                                                         &
    232242        ONLY : rans_mode_parent
    233        
    234    USE random_generator_parallel,                                              &
    235         ONLY:  random_dummy,                                                   &
     243
     244    USE random_generator_parallel,                                             &
     245        ONLY:  init_parallel_random_generator,                                 &
     246               random_dummy,                                                   &
    236247               random_number_parallel,                                         &
    237                random_seed_parallel,                                           &
    238                seq_random_array
     248               random_seed_parallel
    239249
    240250    USE transpose_indices,                                                     &
    … …  
    242252              nzt_x
    243253
     254    USE surface_mod,                                                           &
     255        ONLY:  surf_def_h,                                                     &
     256               surf_lsm_h,                                                     &
     257               surf_usm_h
    244258
    245259    IMPLICIT NONE
    … …  
    253267    LOGICAL ::  parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized
    254268    LOGICAL ::  use_syn_turb_gen = .FALSE.              !< switch to use synthetic turbulence generator
     269    LOGICAL ::  compute_velocity_seeds_local = .TRUE.   !< switch to decide whether velocity seeds are computed locally or if computation
     270                                                        !< is distributed over several processes
    255271
    256272    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator
    … …  
    258274    INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
    259275    INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
    260     INTEGER(iwp) ::  merg               !< maximum length scale (in gp)
    261     INTEGER(iwp) ::  mergp              !< merg + nbgp
     276    INTEGER(iwp) ::  mergp              !< maximum length scale (in gp)
    262277    INTEGER(iwp) ::  nzb_x_stg          !< lower bound of z coordinate (required for transposing z on PEs along x)
    263278    INTEGER(iwp) ::  nzt_x_stg          !< upper bound of z coordinate (required for transposing z on PEs along x)
    … …  
    288303    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwz            !< length scale of w in z direction (in gp)
    289304
     305    INTEGER(isp), DIMENSION(:), ALLOCATABLE   ::  id_rand_xz     !< initial random IDs at xz inflow boundary
     306    INTEGER(isp), DIMENSION(:), ALLOCATABLE   ::  id_rand_yz     !< initial random IDs at yz inflow boundary
     307    INTEGER(isp), DIMENSION(:,:), ALLOCATABLE ::  seq_rand_xz    !< initial random seeds at xz inflow boundary
     308    INTEGER(isp), DIMENSION(:,:), ALLOCATABLE ::  seq_rand_yz    !< initial random seeds at yz inflow boundary
     309
    290310    REAL(wp) ::  blend                    !< value to create gradually and smooth blending of Reynolds stress and length
    291311                                          !< scales above the boundary layer
    … …  
    295315    REAL(wp) ::  length_scale             !< length scale, default is 8 x minimum grid spacing
    296316    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
    297     REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
     317    REAL(wp) ::  dt_stg_call = 0.0_wp     !< time interval for calling synthetic turbulence generator
    298318    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
    299319    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
    … …  
    438458 SUBROUTINE stg_check_parameters
    439459
    440     IMPLICIT NONE
    441 
    442460    IF ( .NOT. use_syn_turb_gen  .AND.  .NOT. rans_mode  .AND.                 &
    443461          nesting_offline )  THEN
    … …  
    526544 SUBROUTINE stg_header ( io )
    527545
    528 
    529     IMPLICIT NONE
    530 
    531546    INTEGER(iwp), INTENT(IN) ::  io   !< Unit of the output file
    532547
    … …  
    562577!------------------------------------------------------------------------------!
    563578 SUBROUTINE stg_init
    564 
    565     IMPLICIT NONE
    566579
    567580    LOGICAL ::  file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist
    … …  
    599612    CALL MPI_BARRIER( comm2d, ierr )
    600613#endif
    601 
     614!
     615!-- Create mpi-datatypes for exchange in case of non-local but distributed
     616!-- computation of the velocity seeds. This option is useful in
     617!-- case large turbulent length scales are presentm, where the computational
     618!-- effort becomes large and need to be parallelized. For parametrized
     619!-- turbulence the length scales are small and computing the velocity seeds
     620!-- locally is faster (no overhead by communication).
     621    IF ( .NOT. compute_velocity_seeds_local )  THEN
    602622#if defined( __parallel )
    603 !
    604 !-- Determine processor decomposition of z-axis along x- and y-direction
    605     nnz = nz / pdims(1)
    606     nzb_x_stg = 1 + myidx * INT( nnz )
    607     nzt_x_stg = ( myidx + 1 ) * INT( nnz )
    608 
    609     IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )             &
    610        nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
    611 
    612     IF ( nesting_offline   .OR.  ( child_domain  .AND.  rans_mode_parent       &
    613                             .AND.  .NOT.  rans_mode ) )  THEN
    614        nnz = nz / pdims(2)
    615        nzb_y_stg = 1 + myidy * INT( nnz )
    616        nzt_y_stg = ( myidy + 1 ) * INT( nnz )
    617 
    618        IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )          &
    619           nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
    620     ENDIF
    621 
    622 !
    623 !-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
    624 !-- velocity seeds
    625     CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
    626     extent = 1 * realsize
    627 !
    628 !-- Set-up MPI datatyp to involve all cores for turbulence generation at yz
    629 !-- layer
    630 !-- stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
    631     CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1],                 &
    632             [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
    633     CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
    634     CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
    635     CALL MPI_TYPE_FREE( newtype, ierr )
    636 
    637     ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
    638     CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyng-nysg+1],     &
    639             [1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
    640     CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
    641     CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
    642     CALL MPI_TYPE_FREE( newtype, ierr )
    643 
    644     ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
    645     ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
    646 
    647     recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
    648     recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
    649 
    650     DO  j = 1, pdims(1)
    651        displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
    652     ENDDO
    653 !
    654 !-- Set-up MPI datatyp to involve all cores for turbulence generation at xz
    655 !-- layer
    656 !-- stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
    657     IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
    658                            .AND.  .NOT.  rans_mode ) )  THEN
    659        CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxrg-nxlg+1],              &
    660                [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
    661        CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
    662        CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
     623!     
     624!--    Determine processor decomposition of z-axis along x- and y-direction
     625       nnz = nz / pdims(1)
     626       nzb_x_stg = 1 + myidx * INT( nnz )
     627       nzt_x_stg = ( myidx + 1 ) * INT( nnz )
     628       
     629       IF ( MOD( nz , pdims(1) ) /= 0  .AND.  myidx == id_stg_right )          &
     630          nzt_x_stg = nzt_x_stg + myidx * ( nnz - INT( nnz ) )
     631       
     632       IF ( nesting_offline   .OR.  ( child_domain  .AND.  rans_mode_parent    &
     633                               .AND.  .NOT.  rans_mode ) )  THEN
     634          nnz = nz / pdims(2)
     635          nzb_y_stg = 1 + myidy * INT( nnz )
     636          nzt_y_stg = ( myidy + 1 ) * INT( nnz )
     637       
     638          IF ( MOD( nz , pdims(2) ) /= 0  .AND.  myidy == id_stg_north )       &
     639             nzt_y_stg = nzt_y_stg + myidy * ( nnz - INT( nnz ) )
     640       ENDIF
     641       
     642!     
     643!--    Define MPI type used in stg_generate_seed_yz to gather vertical splitted
     644!--    velocity seeds
     645       CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
     646       extent = 1 * realsize
     647!     
     648!--    Set-up MPI datatyp to involve all cores for turbulence generation at yz
     649!--    layer
     650!--    stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
     651       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyn-nys+1],                &
     652               [1,nyn-nys+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     653       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
     654       CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
    663655       CALL MPI_TYPE_FREE( newtype, ierr )
    664 
    665        ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
    666        CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxrg-nxlg+1],  &
    667                [1,nxrg-nxlg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
    668        CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
    669        CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
     656       
     657       ! stg_type_yz_small: yz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
     658       CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x_stg-nzb_x_stg+2,nyn-nys+1],    &
     659               [1,nyn-nys+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     660       CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
     661       CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
    670662       CALL MPI_TYPE_FREE( newtype, ierr )
    671 
     663       
    672664       ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
    673        ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
    674 
    675        recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
    676        recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
    677 
    678        DO  j = 1, pdims(2)
    679           displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
     665       ALLOCATE( recv_count_yz(pdims(1)), displs_yz(pdims(1)) )
     666       
     667       recv_count_yz           = nzt_x_stg-nzb_x_stg + 1
     668       recv_count_yz(pdims(1)) = recv_count_yz(pdims(1)) + 1
     669       
     670       DO  j = 1, pdims(1)
     671          displs_yz(j) = 0 + (nzt_x_stg-nzb_x_stg+1) * (j-1)
    680672       ENDDO
    681 
     673!     
     674!--    Set-up MPI datatyp to involve all cores for turbulence generation at xz
     675!--    layer
     676!--    stg_type_xz: xz-slice with vertical bounds nzb:nzt+1
     677       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
     678                              .AND.  .NOT.  rans_mode ) )  THEN
     679          CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nxr-nxl+1],             &
     680                  [1,nxr-nxl+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     681          CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz, ierr )
     682          CALL MPI_TYPE_COMMIT( stg_type_xz, ierr )
     683          CALL MPI_TYPE_FREE( newtype, ierr )
     684       
     685          ! stg_type_yz_small: xz-slice with vertical bounds nzb_x_stg:nzt_x_stg+1
     686          CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_y_stg-nzb_y_stg+2,nxr-nxl+1], &
     687                  [1,nxr-nxl+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
     688          CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_xz_small, ierr )
     689          CALL MPI_TYPE_COMMIT( stg_type_xz_small, ierr )
     690          CALL MPI_TYPE_FREE( newtype, ierr )
     691       
     692          ! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
     693          ALLOCATE( recv_count_xz(pdims(2)), displs_xz(pdims(2)) )
     694       
     695          recv_count_xz           = nzt_y_stg-nzb_y_stg + 1
     696          recv_count_xz(pdims(2)) = recv_count_xz(pdims(2)) + 1
     697       
     698          DO  j = 1, pdims(2)
     699             displs_xz(j) = 0 + (nzt_y_stg-nzb_y_stg+1) * (j-1)
     700          ENDDO
     701       
     702       ENDIF
    682703    ENDIF
    683704
    … …  
    715736               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
    716737               
    717     ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) )
    718     ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) )
     738    ALLOCATE ( dist_xz(nzb:nzt+1,nxl:nxr,3) )
     739    ALLOCATE ( dist_yz(nzb:nzt+1,nys:nyn,3) )
    719740    dist_xz = 0.0_wp
    720741    dist_yz = 0.0_wp
    … …  
    788809!--    Define length scale for the imposed turbulence, which is defined as
    789810!--    8 times the minimum grid spacing
    790        length_scale = 8.0_wp * MIN( dx, dy, MINVAL( dzw ) )
     811       length_scale = 30.0_wp * MIN( dx, dy, MINVAL( dzw ) ) !8.0_wp * MIN( dx, dy, MINVAL( dzw ) )
    791812!
    792813!--    Define constant to gradually decrease length scales and Reynolds stress
    … …  
    837858!
    838859!-- Define the size of the filter functions and allocate them.
    839     merg = 0
     860    mergp = 0
    840861
    841862    ! arrays must be large enough to cover the largest length scale
    … …  
    844865                ABS(nvx(k)), ABS(nvy(k)), ABS(nvz(k)), &
    845866                ABS(nwx(k)), ABS(nwy(k)), ABS(nwz(k))  )
    846        IF ( j > merg )  merg = j
     867       IF ( j > mergp )  mergp = j
    847868    ENDDO
    848869
    849     merg  = 2 * merg
    850     mergp = merg + nbgp
    851 
    852     ALLOCATE ( bux(-merg:merg,nzb:nzt+1),                                      &
    853                buy(-merg:merg,nzb:nzt+1),                                      &
    854                buz(-merg:merg,nzb:nzt+1),                                      &
    855                bvx(-merg:merg,nzb:nzt+1),                                      &
    856                bvy(-merg:merg,nzb:nzt+1),                                      &
    857                bvz(-merg:merg,nzb:nzt+1),                                      &
    858                bwx(-merg:merg,nzb:nzt+1),                                      &
    859                bwy(-merg:merg,nzb:nzt+1),                                      &
    860                bwz(-merg:merg,nzb:nzt+1)  )
     870!     mergp  = 2 * mergp
     871!     mergp = mergp
     872
     873    ALLOCATE ( bux(-mergp:mergp,nzb:nzt+1),                                      &
     874               buy(-mergp:mergp,nzb:nzt+1),                                      &
     875               buz(-mergp:mergp,nzb:nzt+1),                                      &
     876               bvx(-mergp:mergp,nzb:nzt+1),                                      &
     877               bvy(-mergp:mergp,nzb:nzt+1),                                      &
     878               bvz(-mergp:mergp,nzb:nzt+1),                                      &
     879               bwx(-mergp:mergp,nzb:nzt+1),                                      &
     880               bwy(-mergp:mergp,nzb:nzt+1),                                      &
     881               bwz(-mergp:mergp,nzb:nzt+1)  )
    861882
    862883!
    863884!-- Allocate velocity seeds for turbulence at xz-layer
    864     ALLOCATE ( fu_xz( nzb:nzt+1,nxlg:nxrg), fuo_xz(nzb:nzt+1,nxlg:nxrg),       &
    865                fv_xz( nzb:nzt+1,nxlg:nxrg), fvo_xz(nzb:nzt+1,nxlg:nxrg),       &
    866                fw_xz( nzb:nzt+1,nxlg:nxrg), fwo_xz(nzb:nzt+1,nxlg:nxrg)  )
     885    ALLOCATE ( fu_xz( nzb:nzt+1,nxl:nxr), fuo_xz(nzb:nzt+1,nxl:nxr),       &
     886               fv_xz( nzb:nzt+1,nxl:nxr), fvo_xz(nzb:nzt+1,nxl:nxr),       &
     887               fw_xz( nzb:nzt+1,nxl:nxr), fwo_xz(nzb:nzt+1,nxl:nxr)  )
    867888
    868889!
    869890!-- Allocate velocity seeds for turbulence at yz-layer
    870     ALLOCATE ( fu_yz( nzb:nzt+1,nysg:nyng), fuo_yz(nzb:nzt+1,nysg:nyng),       &
    871                fv_yz( nzb:nzt+1,nysg:nyng), fvo_yz(nzb:nzt+1,nysg:nyng),       &
    872                fw_yz( nzb:nzt+1,nysg:nyng), fwo_yz(nzb:nzt+1,nysg:nyng)  )
     891    ALLOCATE ( fu_yz( nzb:nzt+1,nys:nyn), fuo_yz(nzb:nzt+1,nys:nyn),       &
     892               fv_yz( nzb:nzt+1,nys:nyn), fvo_yz(nzb:nzt+1,nys:nyn),       &
     893               fw_yz( nzb:nzt+1,nys:nyn), fwo_yz(nzb:nzt+1,nys:nyn)  )
    873894
    874895    fu_xz  = 0.0_wp
    … …  
    913934          IF ( myidx == id_stg_right )  i = nxr+1
    914935
    915           DO  j = nysg, nyng
     936          DO  j = nys, nyn
    916937             DO  k = nzb, nzt+1
    917938                IF  ( a11(k) > 10E-8_wp )  THEN
    … …  
    944965          IF ( myidy == id_stg_north )  j = nyn+1
    945966
    946           DO  i = nxlg, nxrg
     967          DO  i = nxl, nxr
    947968             DO  k = nzb, nzt+1
    948969!
    … …  
    10391060#endif 
    10401061    ENDIF
     1062!
     1063!-- Initialize random number generator at xz- and yz-layers. Random numbers
     1064!-- are initialized at each core. In case there is only inflow from the left,
     1065!-- it is sufficient to generate only random numbers for the yz-layer, else
     1066!-- random numbers for the xz-layer are also required.
     1067    ALLOCATE ( id_rand_yz(-mergp+nys:nyn+mergp) )
     1068    ALLOCATE ( seq_rand_yz(5,-mergp+nys:nyn+mergp) )
     1069    id_rand_yz  = 0
     1070    seq_rand_yz = 0
     1071
     1072    CALL init_parallel_random_generator( ny, -mergp+nys, nyn+mergp,            &
     1073                                         id_rand_yz, seq_rand_yz )
     1074
     1075    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
     1076                           .AND.  .NOT.  rans_mode ) )  THEN
     1077       ALLOCATE ( id_rand_xz(-mergp+nxl:nxr+mergp) )
     1078       ALLOCATE ( seq_rand_xz(5,-mergp+nxl:nxr+mergp) )
     1079       id_rand_xz  = 0
     1080       seq_rand_xz = 0
     1081
     1082       CALL init_parallel_random_generator( nx, -mergp+nxl, nxr+mergp,         &
     1083                                            id_rand_xz, seq_rand_xz )
     1084    ENDIF
     1085
    10411086
    10421087
    … …  
    10521097 SUBROUTINE stg_filter_func( nxx, bxx )
    10531098
    1054 
    1055     IMPLICIT NONE
    1056 
    10571099    INTEGER(iwp) :: k         !< loop index
    10581100    INTEGER(iwp) :: n_k       !< length scale nXX in height k
    1059     INTEGER(iwp) :: n_k2      !< n_k * 2
    10601101    INTEGER(iwp) :: nf        !< index for length scales
    10611102
    … …  
    10631104    REAL(wp) :: qsi = 1.0_wp  !< minimization factor
    10641105
    1065     INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1)           !< length scale (in gp)
    1066 
    1067     REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1)  !< filter function
     1106    INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  nxx         !< length scale (in gp)
     1107
     1108    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  bxx  !< filter function
    10681109
    10691110
    … …  
    10741115       n_k    = nxx(k)
    10751116       IF ( n_k /= 0 )  THEN
    1076           n_k2 = n_k * 2
    10771117
    10781118!
    10791119!--       ( Eq.10 )^2
    1080           DO  nf = -n_k2, n_k2
     1120          DO  nf = -n_k, n_k
    10811121             bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
    10821122          ENDDO
    … …  
    10851125!--       ( Eq.9 )
    10861126          bdenom = SQRT( bdenom )
    1087           DO  nf = -n_k2, n_k2
     1127          DO  nf = -n_k, n_k
    10881128             bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
    10891129          ENDDO
    … …  
    11011141 SUBROUTINE stg_parin
    11021142
    1103 
    1104     IMPLICIT NONE
    1105 
    11061143    CHARACTER (LEN=80) ::  line   !< dummy string that contains the current line of the parameter file
    11071144
    11081145
    1109     NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
     1146    NAMELIST /stg_par/  dt_stg_adjust,                                         &
     1147                        dt_stg_call,                                           &
     1148                        use_syn_turb_gen,                                      &
     1149                        compute_velocity_seeds_local
    11101150
    11111151    line = ' '
    1112 
    11131152!
    11141153!-- Try to find stg package
    … …  
    11461185 SUBROUTINE stg_rrd_global( found )
    11471186
    1148 
    1149     IMPLICIT NONE
    1150 
    11511187    LOGICAL, INTENT(OUT)  ::  found !< flag indicating if variable was found
    1152 
    11531188
    11541189    found = .TRUE.
    … …  
    11831218 SUBROUTINE stg_wrd_global
    11841219
    1185 
    1186     IMPLICIT NONE
    1187    
    11881220    CALL wrd_write_string( 'time_stg_adjust' )
    11891221    WRITE ( 14 )  time_stg_adjust
    … …  
    12081240!------------------------------------------------------------------------------!
    12091241 SUBROUTINE stg_main
    1210 
    1211     IMPLICIT NONE
    12121242
    12131243    INTEGER(iwp) :: i           !< grid index in x-direction
    … …  
    12521282                  ( child_domain .AND.  rans_mode_parent                       &
    12531283                                 .AND.  .NOT.  rans_mode ) ) )  THEN
    1254 
    12551284          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu_yz, id_stg_left )
    12561285          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv_yz, id_stg_left )
    … …  
    13151344       IF ( stg_call )  THEN
    13161345       
    1317           DO  j = nysg, nyng
     1346          DO  j = nys, nyn
    13181347             DO  k = nzb, nzt + 1
    13191348!     
    … …  
    13481377          IF ( myidx == id_stg_left  )  i = nxl
    13491378          IF ( myidx == id_stg_right )  i = nxr+1       
    1350           DO  j = nysg, nyng
     1379          DO  j = nys, nyn
    13511380             DO  k = nzb+1, nzt + 1
    13521381!     
    … …  
    13611390          IF ( myidx == id_stg_left  )  i = nxl-1
    13621391          IF ( myidx == id_stg_right )  i = nxr+1
    1363           DO  j = nysg, nyng
     1392          DO  j = nys, nyn
    13641393             DO  k = nzb+1, nzt + 1
    13651394!     
    … …  
    14271456          IF ( myidx == id_stg_right )  i = nxr+1                             
    14281457                                                                               
    1429           dist_yz(:,:,1) = ( dist_yz(:,:,1) - mc_factor(1) )                   &
     1458          dist_yz(:,nys:nyn,1) = ( dist_yz(:,nys:nyn,1) - mc_factor(1) )                   &
    14301459                        * MERGE( 1.0_wp, 0.0_wp,                               &
    1431                           BTEST( wall_flags_total_0(:,:,i), 1 ) )             
     1460                          BTEST( wall_flags_total_0(:,nys:nyn,i), 1 ) )             
    14321461                                                                               
    14331462                                                                               
    … …  
    14351464          IF ( myidx == id_stg_right )  i = nxr+1                             
    14361465                                                                               
    1437           dist_yz(:,:,2) = ( dist_yz(:,:,2) - mc_factor(2) )                   &
     1466          dist_yz(:,nys:nyn,2) = ( dist_yz(:,nys:nyn,2) - mc_factor(2) )                   &
    14381467                        * MERGE( 1.0_wp, 0.0_wp,                               &
    1439                           BTEST( wall_flags_total_0(:,:,i), 2 ) )             
     1468                          BTEST( wall_flags_total_0(:,nys:nyn,i), 2 ) )             
    14401469                                                                               
    1441           dist_yz(:,:,3) = ( dist_yz(:,:,3) - mc_factor(3) )                   &
     1470          dist_yz(:,nys:nyn,3) = ( dist_yz(:,nys:nyn,3) - mc_factor(3) )                   &
    14421471                        * MERGE( 1.0_wp, 0.0_wp,                               &
    1443                           BTEST( wall_flags_total_0(:,:,i), 3 ) )
     1472                          BTEST( wall_flags_total_0(:,nys:nyn,i), 3 ) )
    14441473!     
    14451474!--       Add disturbances
    … …  
    14521481!           
    14531482!--             Add disturbance at the inflow
    1454                 DO  j = nysg, nyng
     1483                DO  j = nys, nyn
    14551484                   DO  k = nzb, nzt+1
    14561485                      u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +         &
    … …  
    15101539!--    imposed
    15111540       IF ( stg_call )  THEN
    1512           DO  i = nxlg, nxrg
     1541          DO  i = nxl, nxr
    15131542             DO  k = nzb, nzt + 1
    15141543!         
    … …  
    15441573          IF ( myidy == id_stg_south  ) j = nys
    15451574          IF ( myidy == id_stg_north )  j = nyn+1
    1546           DO  i = nxlg, nxrg
     1575          DO  i = nxl, nxr
    15471576             DO  k = nzb+1, nzt + 1
    15481577!         
    … …  
    15611590          IF ( myidy == id_stg_south  ) j = nys-1
    15621591          IF ( myidy == id_stg_north )  j = nyn+1
    1563           DO  i = nxlg, nxrg
     1592          DO  i = nxl, nxr
    15641593             DO  k = nzb+1, nzt + 1
    15651594!         
    … …  
    16211650          IF ( myidy == id_stg_north )  j = nyn+1
    16221651         
    1623           dist_xz(:,:,2)   = ( dist_xz(:,:,2) - mc_factor(2) )                 &
     1652          dist_xz(:,nxl:nxr,2)   = ( dist_xz(:,nxl:nxr,2) - mc_factor(2) )                 &
    16241653                           * MERGE( 1.0_wp, 0.0_wp,                            &
    1625                              BTEST( wall_flags_total_0(:,j,:), 2 ) )         
     1654                             BTEST( wall_flags_total_0(:,j,nxl:nxr), 2 ) )         
    16261655                                                                               
    16271656                                                                               
    … …  
    16291658          IF ( myidy == id_stg_north )  j = nyn+1                             
    16301659                                                                               
    1631           dist_xz(:,:,1)   = ( dist_xz(:,:,1) - mc_factor(1) )                 &
     1660          dist_xz(:,nxl:nxr,1)   = ( dist_xz(:,nxl:nxr,1) - mc_factor(1) )                 &
    16321661                           * MERGE( 1.0_wp, 0.0_wp,                            &
    1633                              BTEST( wall_flags_total_0(:,j,:), 1 ) )         
     1662                             BTEST( wall_flags_total_0(:,j,nxl:nxr), 1 ) )         
    16341663                                                                               
    1635           dist_xz(:,:,3)   = ( dist_xz(:,:,3) - mc_factor(3) )                 &
     1664          dist_xz(:,nxl:nxr,3)   = ( dist_xz(:,nxl:nxr,3) - mc_factor(3) )                 &
    16361665                           * MERGE( 1.0_wp, 0.0_wp,                            &
    1637                              BTEST( wall_flags_total_0(:,j,:), 3 ) )     
     1666                             BTEST( wall_flags_total_0(:,j,nxl:nxr), 3 ) )     
    16381667!         
    16391668!--       Add disturbances
    … …  
    16731702    ENDIF
    16741703!
     1704!-- Exchange ghost points.
     1705    CALL exchange_horiz( u, nbgp )
     1706    CALL exchange_horiz( v, nbgp )
     1707    CALL exchange_horiz( w, nbgp )
     1708!
    16751709!-- Finally, set time counter for calling STG to zero
    16761710    IF ( stg_call )  time_stg_call = 0.0_wp
    … …  
    16901724!------------------------------------------------------------------------------!
    16911725 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id_left, id_right )
    1692  
    1693  USE pegrid
    1694 
    1695     IMPLICIT NONE
    1696 
    1697     INTEGER(iwp)           :: i           !< grid index x-direction
     1726
    16981727    INTEGER(iwp)           :: id_left     !< core ids at respective boundaries
    16991728    INTEGER(iwp), OPTIONAL :: id_right    !< core ids at respective boundaries
    … …  
    17061735    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y    !< length scale in y-direction
    17071736    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
    1708     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2   !< n_y*2
    1709     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
    17101737
    17111738    REAL(wp) :: nyz_inv         !< inverse of number of grid points in yz-slice
    … …  
    17131740    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
    17141741
    1715     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
    1716     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
    1717    
    1718     REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
    1719     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)             :: f_n     !<  velocity seed
     1742    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1)    :: b_y     !< filter function in y-direction
     1743    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1)    :: b_z     !< filter function in z-direction
     1744   
     1745    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nys:nyn) :: f_n_l   !<  local velocity seed
     1746    REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn)             :: f_n     !<  velocity seed
    17201747   
    17211748    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
    1722     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it_l !< local array of random numbers
    1723 
    17241749!
    17251750!-- Generate random numbers using the parallel random generator.
    … …  
    17321757!-- left boundary, while the right boundary uses the same random numbers
    17331758!-- and thus also computes the same correlation matrix.
    1734     ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
    1735     rand_it   = 0.0_wp
     1759    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nys:nyn+mergp) )
     1760    rand_it = 0.0_wp
    17361761
    17371762    rand_av        = 0.0_wp
    17381763    rand_sigma_inv = 0.0_wp
    1739     nyz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
    1740 !
    1741 !-- Compute and normalize random numbers only on left-boundary ranks.
    1742     IF ( myidx == id_stg_left )  THEN
    1743 !
    1744 !--    Allocate array for local set of random numbers
    1745        ALLOCATE( rand_it_l(nzb:nzt+1,nys:nyn) )
    1746        rand_it_l = 0.0_wp
    1747 
    1748        i = nxl
    1749        DO  j = nys, nyn
    1750 !
    1751 !--       Put the random seeds at grid point j,i
    1752           CALL random_seed_parallel( put=seq_random_array(:, j, i) )
    1753           DO  k = nzb, nzt+1
    1754              CALL random_number_parallel( random_dummy )
    1755              rand_it_l(k,j) = random_dummy
    1756              rand_av        = rand_av + rand_it_l(k,j)
    1757           ENDDO
    1758 !
    1759 !--       Get the new random seeds from last call at grid point j,i
    1760           CALL random_seed_parallel( get=seq_random_array(:, j, i) )
     1764    nyz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 )  &
     1765                                  * ( ny + mergp - ( 0 - mergp ) + 1 ),        &
     1766                                    KIND=wp )
     1767!
     1768!-- Compute and normalize random numbers.
     1769    DO  j = nys - mergp, nyn + mergp
     1770!
     1771!--    Put the random seeds at grid point j
     1772       CALL random_seed_parallel( put=seq_rand_yz(:,j) )
     1773       DO  k = nzb - mergp, nzt + 1 + mergp
     1774          CALL random_number_parallel( random_dummy )
     1775          rand_it(k,j) = random_dummy
    17611776       ENDDO
     1777!
     1778!--    Get the new random seeds from last call at grid point j
     1779       CALL random_seed_parallel( get=seq_rand_yz(:,j) )
     1780    ENDDO
     1781!
     1782!-- For normalization to zero mean, sum-up the global random numers.
     1783!-- To normalize the global set of random numbers,
     1784!-- the inner ghost layers mergp must not be summed-up, else
     1785!-- the random numbers on the ghost layers will be stronger weighted as they
     1786!-- also occur on the inner subdomains.
     1787    DO  j = MERGE( nys, nys - mergp, nys /= 0 ),                              &
     1788            MERGE( nyn, nyn + mergp, nyn /= ny )
     1789       DO  k = nzb - mergp, nzt + 1 + mergp
     1790          rand_av = rand_av + rand_it(k,j)
     1791       ENDDO
     1792    ENDDO
     1793   
    17621794#if defined( __parallel )
    17631795!
    1764 !--    Sum-up the local averages of the random numbers
    1765        CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                 &
    1766                            MPI_SUM, comm1dy, ierr )
     1796!-- Sum-up the local averages of the random numbers
     1797    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                    &
     1798                        MPI_SUM, comm1dy, ierr )
    17671799#endif
    1768        rand_av = rand_av * nyz_inv
    1769 !
    1770 !--    Obtain zero mean
    1771        rand_it_l = rand_it_l - rand_av
    1772 !
    1773 !--    Now, compute the variance
     1800    rand_av = rand_av * nyz_inv
     1801!
     1802!-- Obtain zero mean
     1803    rand_it= rand_it - rand_av
     1804!
     1805!-- Now, compute the variance
     1806    DO  j = MERGE( nys, nys - mergp, nys /= 0 ),                               &
     1807            MERGE( nyn, nyn + mergp, nyn /= ny )
     1808       DO  k = nzb - mergp, nzt + 1 + mergp
     1809          rand_sigma_inv = rand_sigma_inv + rand_it(k,j)**2
     1810       ENDDO
     1811    ENDDO
     1812
     1813#if defined( __parallel )
     1814!
     1815!-- Sum-up the local quadratic averages of the random numbers
     1816    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
     1817                        MPI_SUM, comm1dy, ierr )
     1818#endif
     1819!
     1820!-- Compute standard deviation
     1821    IF ( rand_sigma_inv /= 0.0_wp )  THEN
     1822       rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nyz_inv )
     1823    ELSE
     1824       rand_sigma_inv = 1.0_wp
     1825    ENDIF
     1826!
     1827!-- Normalize with standard deviation to obtain unit variance
     1828    rand_it = rand_it * rand_sigma_inv
     1829
     1830    CALL cpu_log( log_point_s(31), 'STG f_n factors', 'start' )
     1831!
     1832!-- Generate velocity seed following Eq.6 of Xie and Castro (2008). There
     1833!-- are two options. In the first one, the computation of the seeds is
     1834!-- distributed to all processes along the communicator comm1dy and
     1835!-- gathered on the leftmost and, if necessary, on the rightmost process.
     1836!-- For huge length scales the computational effort can become quite huge
     1837!-- (it scales with the turbulent length scales), so that gain by parallelization
     1838!-- exceeds the costs by the subsequent communication.
     1839!-- In the second option, which performs better when the turbulent length scales
     1840!-- are parametrized and thus the loops are smaller, the seeds are computed
     1841!-- locally and no communication is necessary.
     1842    IF ( compute_velocity_seeds_local )  THEN
     1843
     1844       f_n  = 0.0_wp
    17741845       DO  j = nys, nyn
    17751846          DO  k = nzb, nzt+1
    1776              rand_sigma_inv = rand_sigma_inv + rand_it_l(k,j)**2
    1777           ENDDO
    1778        ENDDO
    1779 
    1780 #if defined( __parallel )
    1781 !
    1782 !--    Sum-up the local quadratic averages of the random numbers
    1783        CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
    1784                            MPI_SUM, comm1dy, ierr )
    1785 #endif
    1786 !
    1787 !--    Compute standard deviation
    1788        IF ( rand_sigma_inv /= 0.0_wp )  THEN
    1789           rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nyz_inv )
    1790        ELSE
    1791           rand_sigma_inv = 1.0_wp
    1792        ENDIF
    1793 !
    1794 !--    Normalize with standard deviation to obtain unit variance
    1795        rand_it_l = rand_it_l * rand_sigma_inv
    1796 !
    1797 !--    Copy local random numbers on the global array
    1798        rand_it(nzb:nzt+1,nys:nyn) = rand_it_l(nzb:nzt+1,nys:nyn)
    1799 !
    1800 !--    Deallocate local array
    1801        DEALLOCATE( rand_it_l )
    1802 !
    1803 !--    Now, distribute the final set of random numbers to all mpi ranks located
    1804 !--    on the left boundary. Here, an allreduce with sum reduction is sufficient,
    1805 !--    or, in the non-parallel case, nothing need to be done at all.
    1806 #if defined( __parallel )
    1807        CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,   &
    1808                            MPI_SUM, comm1dy, ierr )
    1809 #endif
    1810     ENDIF
    1811 !
    1812 !-- Finally, distribute the set of random numbers (defined on the leftmost-
    1813 !-- located mpi ranks) to all other mpi ranks. Here, a allreduce with sum
    1814 !-- option is sufficient, because rand_it is zero on all other mpi_ranks.
    1815 !-- Note, the reduce operation is only performed with communicator comm1dx,
    1816 !-- where only 1 rank within the communicator has non-zero random numbers.
    1817 #if defined( __parallel )
    1818     CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,      &
    1819                         MPI_SUM, comm1dx, ierr )
    1820 #endif
    1821 !
    1822 !-- Periodic fill of random numbers in space
    1823     DO  j = 0, ny
    1824        DO  k = 1, mergp
    1825           rand_it(nzb  -k,j) = rand_it(nzt+2-k,j)    ! bottom margin
    1826           rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j)    ! top margin
    1827        ENDDO
    1828     ENDDO
    1829     DO  j = 1, mergp
    1830        DO  k = nzb-mergp, nzt+1+mergp
    1831           rand_it(k,  -j) = rand_it(k,ny-j+1)        ! south margin
    1832           rand_it(k,ny+j) = rand_it(k,   j-1)        ! north margin
    1833        ENDDO
    1834     ENDDO
    1835 
    1836 !
    1837 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
    1838     n_y2 = n_y * 2
    1839     n_z2 = n_z * 2
    1840     f_n_l  = 0.0_wp
    1841 
    1842     DO  j = nysg, nyng
    1843        DO  k = nzb_x_stg, nzt_x_stg+1
    1844           DO  jj = -n_y2(k), n_y2(k)
    1845              DO  kk = -n_z2(k), n_z2(k)
    1846                 f_n_l(k,j) = f_n_l(k,j)                                        &
    1847                            + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
     1847             DO  jj = -n_y(k), n_y(k)
     1848                DO  kk = -n_z(k), n_z(k)
     1849                   f_n(k,j) = f_n(k,j) + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
     1850                ENDDO
    18481851             ENDDO
    18491852          ENDDO
    18501853       ENDDO
    1851     ENDDO
     1854
     1855    ELSE
     1856
     1857       f_n_l  = 0.0_wp
     1858       DO  j = nys, nyn
     1859          DO  k = nzb_x_stg, nzt_x_stg+1
     1860             DO  jj = -n_y(k), n_y(k)
     1861                DO  kk = -n_z(k), n_z(k)
     1862                   f_n_l(k,j) = f_n_l(k,j) + b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
     1863                ENDDO
     1864             ENDDO
     1865          ENDDO
     1866       ENDDO
     1867!
     1868!--    Gather velocity seeds of full subdomain
     1869       send_count = nzt_x_stg - nzb_x_stg + 1
     1870       IF ( nzt_x_stg == nzt )  send_count = send_count + 1
     1871
     1872#if defined( __parallel )
     1873!
     1874!--    Gather the velocity seed matrix on left boundary mpi ranks.
     1875       CALL MPI_GATHERV( f_n_l(nzb_x_stg,nys), send_count, stg_type_yz_small,  &
     1876                         f_n(nzb+1,nys), recv_count_yz, displs_yz, stg_type_yz,&
     1877                         id_left, comm1dx, ierr )
     1878!
     1879!--    If required, gather the same velocity seed matrix on right boundary
     1880!--    mpi ranks (in offline nesting for example).
     1881       IF ( PRESENT( id_right ) )  THEN
     1882          CALL MPI_GATHERV( f_n_l(nzb_x_stg,nys), send_count, stg_type_yz_small,  &
     1883                            f_n(nzb+1,nys), recv_count_yz, displs_yz, stg_type_yz,&
     1884                            id_right, comm1dx, ierr )
     1885       ENDIF
     1886#else
     1887       f_n(nzb+1:nzt+1,nys:nyn) = f_n_l(nzb_x_stg:nzt_x_stg+1,nys:nyn)
     1888#endif
     1889
     1890    ENDIF
    18521891
    18531892    DEALLOCATE( rand_it )
    1854 !
    1855 !-- Gather velocity seeds of full subdomain
    1856     send_count = nzt_x_stg - nzb_x_stg + 1
    1857     IF ( nzt_x_stg == nzt )  send_count = send_count + 1
    1858 
    1859 #if defined( __parallel )
    1860 !
    1861 !-- Gather the velocity seed matrix on left boundary mpi ranks.
    1862     CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,     &
    1863                       f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,   &
    1864                       id_left, comm1dx, ierr )
    1865 !
    1866 !-- If required, gather the same velocity seed matrix on right boundary
    1867 !-- mpi ranks (in offline nesting for example).
    1868     IF ( PRESENT( id_right ) )  THEN
    1869        CALL MPI_GATHERV( f_n_l(nzb_x_stg,nysg), send_count, stg_type_yz_small,  &
    1870                          f_n(nzb+1,nysg), recv_count_yz, displs_yz, stg_type_yz,&
    1871                          id_right, comm1dx, ierr )
    1872     ENDIF
    1873 #else
    1874     f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x_stg:nzt_x_stg+1,nysg:nyng)
    1875 #endif
    1876 
     1893
     1894    CALL cpu_log( log_point_s(31), 'STG f_n factors', 'stop' )
    18771895
    18781896 END SUBROUTINE stg_generate_seed_yz
    … …  
    18901908 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north )
    18911909
    1892     IMPLICIT NONE
    1893 
    18941910    INTEGER(iwp) :: i           !< loop index in x-direction
    18951911    INTEGER(iwp) :: id_north    !< core ids at respective boundaries
    18961912    INTEGER(iwp) :: id_south    !< core ids at respective boundaries
    18971913    INTEGER(iwp) :: ii          !< loop index in x-direction
    1898     INTEGER(iwp) :: j           !< grid index y-direction
    18991914    INTEGER(iwp) :: k           !< loop index in z-direction
    19001915    INTEGER(iwp) :: kk          !< loop index in z-direction
    … …  
    19031918    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
    19041919    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
    1905     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_x*2
    1906     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
    19071920
    19081921    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
    … …  
    19101923    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
    19111924
    1912     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
    1913     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
    1914    
    1915     REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
    1916     REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)             :: f_n     !<  velocity seed
     1925    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1)    :: b_x     !< filter function in x-direction
     1926    REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1)    :: b_z     !< filter function in z-direction
     1927   
     1928    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) :: f_n_l   !<  local velocity seed
     1929    REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr)             :: f_n     !<  velocity seed
    19171930
    19181931    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
    1919     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it_l !< local array of random numbers
    19201932
    19211933!
    … …  
    19291941!-- left boundary, while the right boundary uses the same random numbers
    19301942!-- and thus also computes the same correlation matrix.
    1931     ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:nx+mergp) )
    1932     rand_it   = 0.0_wp
     1943    ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nxl:nxr+mergp) )
     1944    rand_it = 0.0_wp
    19331945
    19341946    rand_av        = 0.0_wp
    19351947    rand_sigma_inv = 0.0_wp
    1936     nxz_inv        = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( nx+1 ), KIND=wp )
    1937 !
    1938 !-- Compute and normalize random numbers only on south-boundary ranks.
    1939     IF ( myidy == id_stg_south )  THEN
    1940 !
    1941 !--    Allocate array for local set of random numbers
    1942        ALLOCATE( rand_it_l(nzb:nzt+1,nxl:nxr) )
    1943        rand_it_l = 0.0_wp
    1944 
    1945        j = nys
    1946        DO  i = nxl, nxr
    1947 !
    1948 !--       Put the random seeds at grid point j,i
    1949           CALL random_seed_parallel( put=seq_random_array(:, j, i) )
    1950           DO  k = nzb, nzt+1
    1951              CALL random_number_parallel( random_dummy )
    1952              rand_it_l(k,i) = random_dummy
    1953              rand_av = rand_av + rand_it_l(k,i)
    1954           ENDDO
    1955 !
    1956 !--       Get the new random seeds from last call at grid point j,i
    1957           CALL random_seed_parallel( get=seq_random_array(:, j, i) )
     1948    nxz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 )  &
     1949                                  * ( nx + mergp - ( 0 - mergp ) +1 ),         &
     1950                                    KIND=wp )
     1951!
     1952!-- Compute and normalize random numbers.
     1953    DO  i = nxl - mergp, nxr + mergp
     1954!
     1955!--    Put the random seeds at grid point ii
     1956       CALL random_seed_parallel( put=seq_rand_xz(:,i) )
     1957       DO  k = nzb - mergp, nzt + 1 + mergp
     1958          CALL random_number_parallel( random_dummy )
     1959          rand_it(k,i) = random_dummy
    19581960       ENDDO
     1961!
     1962!--    Get the new random seeds from last call at grid point ii
     1963       CALL random_seed_parallel( get=seq_rand_xz(:,i) )
     1964    ENDDO
     1965!
     1966!-- For normalization to zero mean, sum-up the global random numers.
     1967!-- To normalize the global set of random numbers,
     1968!-- the inner ghost layers mergp must not be summed-up, else
     1969!-- the random numbers on the ghost layers will be stronger weighted as they
     1970!-- also occur on the inner subdomains.
     1971    DO  i = MERGE( nxl, nxl - mergp, nxl /= 0 ),                              &
     1972            MERGE( nxr, nxr + mergp, nxr /= nx )
     1973       DO  k = nzb - mergp, nzt + 1 + mergp
     1974          rand_av = rand_av + rand_it(k,i)
     1975       ENDDO
     1976    ENDDO
     1977   
    19591978#if defined( __parallel )
    1960        CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                 &
    1961                            MPI_SUM, comm1dx, ierr )
     1979!
     1980!-- Sum-up the local averages of the random numbers
     1981    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_av, 1, MPI_REAL,                    &
     1982                        MPI_SUM, comm1dx, ierr )
    19621983#endif
    1963        rand_av = rand_av * nxz_inv
    1964 !
    1965 !--    Obtain zero mean
    1966        rand_it_l = rand_it_l - rand_av
    1967 !
    1968 !--    Now, compute the variance
     1984    rand_av = rand_av * nxz_inv
     1985!
     1986!-- Obtain zero mean
     1987    rand_it= rand_it - rand_av
     1988!
     1989!-- Now, compute the variance
     1990    DO  i = MERGE( nxl, nxl - mergp, nxl /= 0 ),                               &
     1991            MERGE( nxr, nxr + mergp, nxr /= nx )
     1992       DO  k = nzb - mergp, nzt + 1 + mergp
     1993          rand_sigma_inv = rand_sigma_inv + rand_it(k,i)**2
     1994       ENDDO
     1995    ENDDO
     1996
     1997#if defined( __parallel )
     1998!
     1999!-- Sum-up the local quadratic averages of the random numbers
     2000    CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
     2001                        MPI_SUM, comm1dx, ierr )
     2002#endif
     2003!
     2004!-- Compute standard deviation
     2005    IF ( rand_sigma_inv /= 0.0_wp )  THEN
     2006       rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nxz_inv )
     2007    ELSE
     2008       rand_sigma_inv = 1.0_wp
     2009    ENDIF
     2010!
     2011!-- Normalize with standard deviation to obtain unit variance
     2012    rand_it = rand_it * rand_sigma_inv
     2013
     2014    CALL cpu_log( log_point_s(31), 'STG f_n factors', 'start' )
     2015!
     2016!-- Generate velocity seed following Eq.6 of Xie and Castro (2008). There
     2017!-- are two options. In the first one, the computation of the seeds is
     2018!-- distributed to all processes along the communicator comm1dx and
     2019!-- gathered on the southmost and, if necessary, on the northmost process.
     2020!-- For huge length scales the computational effort can become quite huge
     2021!-- (it scales with the turbulent length scales), so that gain by parallelization
     2022!-- exceeds the costs by the subsequent communication.
     2023!-- In the second option, which performs better when the turbulent length scales
     2024!-- are parametrized and thus the loops are smaller, the seeds are computed
     2025!-- locally and no communication is necessary.
     2026    IF ( compute_velocity_seeds_local )  THEN
     2027
     2028       f_n  = 0.0_wp
    19692029       DO  i = nxl, nxr
    19702030          DO  k = nzb, nzt+1
    1971              rand_sigma_inv = rand_sigma_inv + rand_it_l(k,i)**2
    1972           ENDDO
    1973        ENDDO
    1974 
    1975 #if defined( __parallel )
    1976        CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_sigma_inv, 1, MPI_REAL,          &
    1977                            MPI_SUM, comm1dx, ierr )
    1978 #endif
    1979 !
    1980 !--    Compute standard deviation
    1981        IF ( rand_sigma_inv /= 0.0_wp )  THEN
    1982           rand_sigma_inv = 1.0_wp / SQRT( rand_sigma_inv * nxz_inv )
    1983        ELSE
    1984           rand_sigma_inv = 1.0_wp
    1985        ENDIF
    1986 !
    1987 !--    Normalize with standard deviation to obtain unit variance
    1988        rand_it_l = rand_it_l * rand_sigma_inv
    1989 !
    1990 !--    Copy local random numbers on the global array
    1991        rand_it(nzb:nzt+1,nxl:nxr) = rand_it_l(nzb:nzt+1,nxl:nxr)
    1992 !
    1993 !--    Deallocate local array
    1994        DEALLOCATE( rand_it_l )
    1995 !
    1996 !--    Now, distribute the final set of random numbers to all mpi ranks located
    1997 !--    on the south boundary. Here, an allreduce with sum reduction is sufficient.
    1998 #if defined( __parallel )
    1999        CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,   &
    2000                            MPI_SUM, comm1dx, ierr )
    2001 #endif
    2002     ENDIF
    2003 !
    2004 !-- Finally, distribute the set of random numbers (defined on the 
    2005 !-- southmost-located mpi ranks) to all other mpi ranks. Here, a allreduce
    2006 !-- with sum option is sufficient, because rand_it is zero on all other
    2007 !-- mpi_ranks. Note, the reduce operation is only performed with communicator
    2008 !-- comm1dy, where only 1 rank within the communicator has non zero random numbers.
    2009 #if defined( __parallel )
    2010     CALL MPI_ALLREDUCE( MPI_IN_PLACE, rand_it, SIZE( rand_it ), MPI_REAL,      &
    2011                         MPI_SUM, comm1dy, ierr )
    2012 #endif
    2013 !
    2014 !-- Periodic fill of random number in space
    2015     DO  i = 0, nx
    2016        DO  k = 1, mergp
    2017           rand_it(nzb-k,i)   = rand_it(nzt+2-k,i)    ! bottom margin
    2018           rand_it(nzt+1+k,i) = rand_it(nzb+k-1,i)    ! top margin
    2019        ENDDO
    2020     ENDDO
    2021     DO  i = 1, mergp
    2022        DO  k = nzb-mergp, nzt+1+mergp
    2023           rand_it(k,-i)   = rand_it(k,nx-i+1)        ! left margin
    2024           rand_it(k,nx+i) = rand_it(k,i-1)           ! right margin
    2025        ENDDO
    2026     ENDDO
    2027 !
    2028 !-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
    2029     n_x2 = n_x * 2
    2030     n_z2 = n_z * 2
    2031     f_n_l  = 0.0_wp
    2032 
    2033     DO  i = nxlg, nxrg
    2034        DO  k = nzb_y_stg, nzt_y_stg+1
    2035           DO  ii = -n_x2(k), n_x2(k)
    2036              DO  kk = -n_z2(k), n_z2(k)
    2037                 f_n_l(k,i) = f_n_l(k,i)                                        &
    2038                            + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
     2031             DO  ii = -n_x(k), n_x(k)
     2032                DO  kk = -n_z(k), n_z(k)
     2033                   f_n(k,i) = f_n(k,i) + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
     2034                ENDDO
    20392035             ENDDO
    20402036          ENDDO
    20412037       ENDDO
    2042     ENDDO
     2038
     2039    ELSE
     2040
     2041       f_n_l  = 0.0_wp
     2042       DO  i = nxl, nxr
     2043          DO  k = nzb_y_stg, nzt_y_stg+1
     2044             DO  ii = -n_x(k), n_x(k)
     2045                DO  kk = -n_z(k), n_z(k)
     2046                   f_n_l(k,i) = f_n_l(k,i) + b_x(ii,k) * b_z(kk,k) * rand_it(k+kk,i+ii)
     2047                ENDDO
     2048             ENDDO
     2049          ENDDO
     2050       ENDDO
     2051!
     2052!--    Gather velocity seeds of full subdomain
     2053       send_count = nzt_y_stg - nzb_y_stg + 1
     2054       IF ( nzt_y_stg == nzt )  send_count = send_count + 1
     2055
     2056#if defined( __parallel )
     2057!
     2058!--    Gather the processed velocity seed on south boundary mpi ranks.
     2059       CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxl), send_count, stg_type_xz_small,   &
     2060                         f_n(nzb+1,nxl), recv_count_xz, displs_xz, stg_type_xz, &
     2061                         id_south, comm1dy, ierr )
     2062!
     2063!--    Gather the processed velocity seed on north boundary mpi ranks.
     2064       CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxl), send_count, stg_type_xz_small,   &
     2065                         f_n(nzb+1,nxl), recv_count_xz, displs_xz, stg_type_xz, &
     2066                         id_north, comm1dy, ierr )
     2067#else
     2068       f_n(nzb+1:nzt+1,nxl:nxr) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxl:nxr)
     2069#endif
     2070
     2071    ENDIF
    20432072
    20442073    DEALLOCATE( rand_it )
    20452074
    2046 !
    2047 !-- Gather velocity seeds of full subdomain
    2048     send_count = nzt_y_stg - nzb_y_stg + 1
    2049     IF ( nzt_y_stg == nzt )  send_count = send_count + 1
    2050 
    2051 #if defined( __parallel )
    2052 !
    2053 !-- Gather the processed velocity seed on south boundary mpi ranks.
    2054     CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,   &
    2055                       f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
    2056                       id_south, comm1dy, ierr )
    2057 !
    2058 !-- Gather the processed velocity seed on north boundary mpi ranks.
    2059     CALL MPI_GATHERV( f_n_l(nzb_y_stg,nxlg), send_count, stg_type_xz_small,   &
    2060                       f_n(nzb+1,nxlg), recv_count_xz, displs_xz, stg_type_xz, &
    2061                       id_north, comm1dy, ierr )
    2062 #else
    2063     f_n(nzb+1:nzt+1,nxlg:nxrg) = f_n_l(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg)
    2064 #endif
     2075    CALL cpu_log( log_point_s(31), 'STG f_n factors', 'stop' )
    20652076
    20662077 END SUBROUTINE stg_generate_seed_xz
    … …  
    20762087 SUBROUTINE parametrize_reynolds_stress
    20772088
    2078     IMPLICIT NONE
    2079 
    20802089    INTEGER(iwp) :: k            !< loop index in z-direction
    20812090   
    … …  
    20932102!
    20942103!--    u'u' and v'v'. Assume isotropy. Note, add a small negative number
    2095 !--    to the denominator, else the merge-function can crash if scale_l is
     2104!--    to the denominator, else the mergpe-function can crash if scale_l is
    20962105!--    zero.
    20972106       r11(k) = scale_us**2 * (                                                &
    … …  
    21402149!------------------------------------------------------------------------------!
    21412150 SUBROUTINE calc_coeff_matrix
    2142 
    2143     IMPLICIT NONE
    21442151
    21452152    INTEGER(iwp) :: k   !< loop index in z-direction
    … …  
    21952202 SUBROUTINE stg_adjust
    21962203
    2197     IMPLICIT NONE
    2198 
    2199 
    22002204    IF ( debug_output_timestep )  CALL debug_message( 'stg_adjust', 'start' )
    22012205!
    … …  
    22512255 SUBROUTINE calc_length_and_time_scale
    22522256
    2253     IMPLICIT NONE
    2254 
    2255 
    22562257    INTEGER(iwp) ::  k !< loop index in z-direction
    22572258   
    … …  
    22852286!--    Assume isotropic turbulence length scales
    22862287       nux(k) = MAX( INT( length_scale * ddx     ), 1 ) * blend
    2287        nuy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend       
     2288       nuy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend
    22882289       nvx(k) = MAX( INT( length_scale * ddx     ), 1 ) * blend
    2289        nvy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend       
     2290       nvy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend
    22902291       nwx(k) = MAX( INT( length_scale * ddx     ), 1 ) * blend
    2291        nwy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend       
     2292       nwy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend
    22922293!
    22932294!--    Along the vertical direction limit the length scale further by the
    … …  
    23412342!------------------------------------------------------------------------------!
    23422343 SUBROUTINE calc_scaling_variables
    2343        
    2344     USE surface_mod,                                                           &
    2345         ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
    2346 
    2347     IMPLICIT NONE
    23482344
    23492345    INTEGER(iwp) :: i            !< loop index in x-direction
Note: See TracChangeset for help on using the changeset viewer.
                                                                                                                                                                                                                                                                                                                                                                               
  | Impressum | ©Leibniz Universität Hannover |