Changeset 4022 for palm/trunk


Ignore:
Timestamp:
Jun 12, 2019 11:52:39 AM (5 years ago)
Author:
suehring
Message:

Synthetic turbulence generator: Revise bias correction of imposed perturbations (correction via volume flow can create instabilities in case the mean volume flow is close to zero); Introduce lower limits in calculation of coefficient matrix, else the calculation may become numerically unstable; Impose perturbations every timestep, even though no new set of perturbations is generated in case dt_stg_call /= dt_3d; Implement a gradual decrease of Reynolds stress and length scales above ABL height (within 1 length scale above ABL depth to 1/10) rather than an abrupt decline; Bugfix in non-nested case: use ABL height for parametrized turbulence; Offline nesting: Rename subroutine for ABL height calculation and make it public; Change top boundary condition for pressure back to Neumann

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

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

    r3987 r4022  
    2525! -----------------
    2626! $Id$
     27! Detection of boundary-layer depth in stable boundary layer on basis of
     28! boundary data improved
     29! Routine for boundary-layer depth calculation renamed and made public
     30!
     31! 3987 2019-05-22 09:52:13Z kanani
    2732! Introduce alternative switch for debug output during timestepping
    2833!
     
    117122!
    118123!-- Public subroutines
    119     PUBLIC nesting_offl_bc, nesting_offl_check_parameters, nesting_offl_header,&
    120            nesting_offl_init, nesting_offl_mass_conservation, nesting_offl_parin
     124    PUBLIC nesting_offl_bc,                                                    &
     125           nesting_offl_calc_zi,                                               &
     126           nesting_offl_check_parameters,                                      &
     127           nesting_offl_header,                                                &
     128           nesting_offl_init,                                                  &
     129           nesting_offl_mass_conservation,                                     &
     130           nesting_offl_parin
    121131!
    122132!-- Public variables
     
    126136       MODULE PROCEDURE nesting_offl_bc
    127137    END INTERFACE nesting_offl_bc
     138   
     139    INTERFACE nesting_offl_calc_zi
     140       MODULE PROCEDURE nesting_offl_calc_zi
     141    END INTERFACE nesting_offl_calc_zi
    128142   
    129143    INTERFACE nesting_offl_check_parameters
     
    806820!--    Further, adjust Rayleigh damping height in case of time-changing conditions.
    807821!--    Therefore, calculate boundary-layer depth first.
    808        CALL calc_zi
     822       CALL nesting_offl_calc_zi
    809823       CALL adjust_sponge_layer
    810824   
     
    833847!> bulk-Richardson criterion.
    834848!------------------------------------------------------------------------------!
    835     SUBROUTINE calc_zi
     849    SUBROUTINE nesting_offl_calc_zi
    836850       
    837851       USE basic_constants_and_equations_mod,                                  &
     
    845859       IMPLICIT NONE
    846860
    847        INTEGER(iwp) :: i            !< loop index in x-direction
    848        INTEGER(iwp) :: j            !< loop index in y-direction
    849        INTEGER(iwp) :: k            !< loop index in z-direction
    850        INTEGER(iwp) :: k_surface    !< topography top index in z-direction
     861       INTEGER(iwp) :: i                            !< loop index in x-direction
     862       INTEGER(iwp) :: j                            !< loop index in y-direction
     863       INTEGER(iwp) :: k                            !< loop index in z-direction
     864       INTEGER(iwp) :: k_max_loc                    !< index of maximum wind speed along z-direction
     865       INTEGER(iwp) :: k_surface                    !< topography top index in z-direction
     866       INTEGER(iwp) :: num_boundary_gp_non_cyclic   !< number of non-cyclic boundaries, used for averaging ABL depth
     867       INTEGER(iwp) :: num_boundary_gp_non_cyclic_l !< number of non-cyclic boundaries, used for averaging ABL depth
    851868   
    852869       REAL(wp) ::  ri_bulk                 !< bulk Richardson number
     
    854871       REAL(wp) ::  topo_max                !< maximum topography level in model domain
    855872       REAL(wp) ::  topo_max_l              !< maximum topography level in subdomain
    856        REAL(wp) ::  u_comp                  !< u-component
    857        REAL(wp) ::  v_comp                  !< v-component
    858873       REAL(wp) ::  vpt_surface             !< near-surface virtual potential temperature
    859874       REAL(wp) ::  zi_l                    !< mean boundary-layer depth on subdomain   
     
    861876   
    862877       REAL(wp), DIMENSION(nzb:nzt+1) ::  vpt_col !< vertical profile of virtual potential temperature at (j,i)-grid point
     878       REAL(wp), DIMENSION(nzb:nzt+1) ::  uv_abs  !< vertical profile of horizontal wind speed at (j,i)-grid point
    863879
    864880       
     
    867883!--    Start with the left and right boundaries.
    868884       zi_l      = 0.0_wp
     885       num_boundary_gp_non_cyclic_l = 0
    869886       IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
     887!
     888!--       Sum-up and store number of boundary grid points used for averaging
     889!--       ABL depth
     890          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
     891                                         nxr - nxl + 1
    870892!
    871893!--       Determine index along x. Please note, index indicates boundary
     
    898920!--          the boundary, where velocity inside the building is zero
    899921!--          (k_surface is the index of the lowest upward-facing surface).
     922             uv_abs(:) = SQRT( MERGE( u(:,j,i+1), u(:,j,i),                   &
     923                                      bc_dirichlet_l )**2 +                   &
     924                               v(:,j,i)**2 )
     925!
     926!--          Determine index of the maximum wind speed                               
     927             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
     928
    900929             zi_local = 0.0_wp
    901930             DO  k = k_surface+1, nzt
    902                 u_comp = MERGE( u(k,j,i+1), u(k,j,i), bc_dirichlet_l )
    903                 v_comp = v(k,j,i)
    904931                ri_bulk = zu(k) * g / vpt_surface *                            &
    905932                          ( vpt_col(k) - vpt_surface ) /                       &
    906                           ( u_comp * u_comp + v_comp * v_comp + 1E-5_wp )
    907                        
    908                 IF ( zi_local == 0.0_wp  .AND.  ri_bulk > ri_bulk_crit )       &
     933                          ( uv_abs(k) + 1E-5_wp )
     934!
     935!--             Check if critical Richardson number is exceeded. Further, check
     936!--             if there is a maxium in the wind profile in order to detect also
     937!--             ABL heights in the stable boundary layer.
     938                IF ( zi_local == 0.0_wp  .AND.                                 &
     939                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
    909940                   zi_local = zu(k)           
    910941             ENDDO
     
    920951!--    Do the same at the north and south boundaries.
    921952       IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
     953       
     954          num_boundary_gp_non_cyclic_l = num_boundary_gp_non_cyclic_l +        &
     955                                         nxr - nxl + 1
    922956       
    923957          j = MERGE( -1, nyn + 1, bc_dirichlet_s )
     
    935969             ENDIF
    936970         
     971             uv_abs(:) = SQRT( u(:,j,i)**2 +                                   &
     972                               MERGE( v(:,j+1,i), v(:,j,i),                    &
     973                               bc_dirichlet_s )**2 )
     974!
     975!--          Determine index of the maximum wind speed
     976             k_max_loc = MAXLOC( uv_abs(:), DIM = 1 ) - 1
     977         
    937978             zi_local = 0.0_wp
    938979             DO  k = k_surface+1, nzt               
    939                 u_comp = u(k,j,i)
    940                 v_comp = MERGE( v(k,j+1,i), v(k,j,i), bc_dirichlet_s )
    941980                ri_bulk = zu(k) * g / vpt_surface *                            &
    942981                       ( vpt_col(k) - vpt_surface ) /                          &
    943                        ( u_comp * u_comp + v_comp * v_comp + 1E-5_wp )
    944                        
    945                 IF ( zi_local == 0.0_wp  .AND.  ri_bulk > 0.25_wp )            &
     982                       ( uv_abs(k) + 1E-5_wp )
     983!
     984!--             Check if critical Richardson number is exceeded. Further, check
     985!--             if there is a maxium in the wind profile in order to detect also
     986!--             ABL heights in the stable boundary layer.                       
     987                IF ( zi_local == 0.0_wp  .AND.                                 &
     988                     ( ri_bulk > ri_bulk_crit  .OR.  k == k_max_loc ) )        &
    946989                   zi_local = zu(k)   
    947990             ENDDO
     
    955998       CALL MPI_ALLREDUCE( zi_l, zi_ribulk, 1, MPI_REAL, MPI_SUM,              &
    956999                           comm2d, ierr )
     1000       CALL MPI_ALLREDUCE( num_boundary_gp_non_cyclic_l,                       &
     1001                           num_boundary_gp_non_cyclic,                         &
     1002                           1, MPI_INTEGER, MPI_SUM, comm2d, ierr )
    9571003#else
    9581004       zi_ribulk = zi_l
     1005       num_boundary_gp_non_cyclic = num_boundary_gp_non_cyclic_l
    9591006#endif
    960        zi_ribulk = zi_ribulk / REAL( 2 * nx + 2 * ny, KIND = wp )
     1007       zi_ribulk = zi_ribulk / REAL( num_boundary_gp_non_cyclic, KIND = wp )
    9611008!
    9621009!--    Finally, check if boundary layer depth is not below the any topography.
     
    9681015       
    9691016#if defined( __parallel )
    970        CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,            &
     1017       CALL MPI_ALLREDUCE( topo_max_l, topo_max, 1, MPI_REAL, MPI_MAX,         &
    9711018                           comm2d, ierr )
    9721019#else
    9731020       topo_max     = topo_max_l
    9741021#endif
    975        zi_ribulk = MAX( zi_ribulk, topo_max )
    976        
    977     END SUBROUTINE calc_zi
     1022!        zi_ribulk = MAX( zi_ribulk, topo_max )
     1023       
     1024    END SUBROUTINE nesting_offl_calc_zi
    9781025   
    9791026   
     
    13191366!--    boundary data. This is requiered for initialize the synthetic turbulence
    13201367!--    generator correctly.
    1321        CALL calc_zi
     1368       CALL nesting_offl_calc_zi
    13221369       
    13231370!
  • palm/trunk/SOURCE/parin.f90

    r4017 r4022  
    2525! -----------------
    2626! $Id$
     27! Change default top boundary condition for pressure to Neumann in offline
     28! nesting case
     29!
     30! 4017 2019-06-06 12:16:46Z schwenkel
    2731! Introduce alternative switch for debug output during timestepping
    2832!
     
    983987             bc_pt_t  = 'nesting_offline'
    984988             bc_q_t   = 'nesting_offline'
    985            !  bc_s_t   = 'nesting_offline'  ! scalar boundary condition is not clear
     989           !  bc_s_t   = 'nesting_offline'  ! scalar boundary condition is not clear yet
    986990           !  bc_cs_t  = 'nesting_offline'  ! same for chemical species
    987 !
    988 !--          For the pressure set Dirichlet conditions, in contrast to the
    989 !--          self nesting. This gives less oscilations within the
    990 !--          free atmosphere since the pressure solver has more degrees of
    991 !--          freedom. In constrast to the self nesting, this might be
    992 !--          justified since the top boundary is far away from the domain
    993 !--          of interest.
    994              bc_p_t   = 'dirichlet' !'neumann'
     991             bc_p_t   = 'neumann'
    995992          ENDIF
    996993
  • palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90

    r3987 r4022  
    2525! -----------------
    2626! $Id$
     27! Several bugfixes and improvements
     28! - revise bias correction of the imposed perturbations (correction via volume
     29!   flow can create instabilities in case the mean volume flow is close to zero)
     30! - introduce lower limits in calculation of coefficient matrix, else the
     31!   calculation may become numerically unstable
     32! - impose perturbations every timestep, even though no new set of perturbations
     33!   is generated in case dt_stg_call /= dt_3d
     34! - Implement a gradual decrease of Reynolds stress and length scales above
     35!   ABL height (within 1 length scale above ABL depth to 1/10) rather than a
     36!   discontinuous decrease
     37! - Bugfix in non-nested case: use ABL height for parametrized turbulence
     38!
     39! 3987 2019-05-22 09:52:13Z kanani
    2740! Introduce alternative switch for debug output during timestepping
    2841!
     
    206219
    207220    USE arrays_3d,                                                             &
    208         ONLY:  mean_inflow_profiles, q, pt, u, v, w, zu, zw
     221        ONLY:  dzw,                                                            &
     222               ddzw,                                                           &
     223               drho_air,                                                       &
     224               mean_inflow_profiles,                                           &
     225               q,                                                              &
     226               pt,                                                             &
     227               u,                                                              &
     228               u_init,                                                         &
     229               v,                                                              &
     230               v_init,                                                         &
     231               w,                                                              &
     232               zu,                                                             &
     233               zw
    209234
    210235    USE basic_constants_and_equations_mod,                                     &
    211         ONLY:  g, kappa, pi
     236        ONLY:  g,                                                              &
     237               kappa,                                                          &
     238               pi
    212239
    213240    USE control_parameters,                                                    &
    214         ONLY:  debug_output_timestep,                                          &
     241        ONLY:  bc_lr,                                                          &
     242               bc_ns,                                                          &
     243               child_domain,                                                   &
     244               coupling_char,                                                  &
     245               debug_output_timestep,                                          &
     246               dt_3d,                                                          &
     247               e_init,                                                         &
    215248               initializing_actions,                                           &
     249               intermediate_timestep_count,                                    &
     250               intermediate_timestep_count_max,                                &
     251               length,                                                         &
    216252               message_string,                                                 &
     253               nesting_offline,                                                &
    217254               num_mean_inflow_profiles,                                       &
    218                syn_turb_gen
     255               rans_mode,                                                      &
     256               restart_string,                                                 &
     257               syn_turb_gen,                                                   &
     258               time_since_reference_point,                                     &
     259               turbulent_inflow
     260               
     261    USE grid_variables,                                                        &
     262        ONLY:  ddx,                                                            &
     263               ddy,                                                            &
     264               dx,                                                             &
     265               dy
    219266
    220267    USE indices,                                                               &
    221         ONLY:  nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg
     268        ONLY:  nbgp,                                                           &
     269               nz,                                                             &
     270               nzb,                                                            &
     271               nzt,                                                            &
     272               nx,                                                             &
     273               nxl,                                                            &
     274               nxlu,                                                           &
     275               nxlg,                                                           &
     276               nxr,                                                            &
     277               nxrg,                                                           &
     278               ny,                                                             &
     279               nys,                                                            &
     280               nysv,                                                           &
     281               nyn,                                                            &
     282               nyng,                                                           &
     283               nysg,                                                           &
     284               wall_flags_0
    222285
    223286    USE kinds
     
    228291
    229292    USE nesting_offl_mod,                                                      &
    230         ONLY:  zi_ribulk
     293        ONLY:  nesting_offl_calc_zi,                                           &
     294               zi_ribulk
    231295
    232296    USE pegrid,                                                                &
    233         ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
     297        ONLY:  comm1dx,                                                        &
     298               comm1dy,                                                        &
     299               comm2d,                                                         &
     300               ierr,                                                           &
     301               myidx,                                                          &
     302               myidy,                                                          &
     303               pdims
     304       
     305    USE pmc_interface,                                                         &
     306        ONLY : rans_mode_parent
    234307
    235308    USE transpose_indices,                                                     &
    236         ONLY: nzb_x, nzt_x
     309        ONLY: nzb_x,                                                           &
     310              nzt_x
    237311
    238312
     
    263337    INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
    264338
     339    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections,
     340                                                      !< required for bias correction of imposed perturbations
     341    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_yz = 0 !< number of non-topography grid points at yz cross-sections,
     342                                                      !< required for bias correction of imposed perturbations
     343   
    265344    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_xz      !< displacement for MPI_GATHERV
    266345    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_xz  !< receive count for MPI_GATHERV
     
    278357    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
    279358
     359    REAL(wp) ::  blend                    !< value to create gradually and smooth blending of Reynolds stress and length
     360                                          !< scales above the boundary layer
     361    REAL(wp) ::  blend_coeff = -2.3_wp    !< coefficient used to ensure that blending functions decreases to 1/10 after
     362                                          !< one length scale above ABL top
     363    REAL(wp) ::  d_l                      !< blend_coeff/length_scale
     364    REAL(wp) ::  length_scale             !< length scale, default is 8 x minimum grid spacing
    280365    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
    281366    REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
    282     REAL(wp) ::  mc_factor = 1.0_wp       !< mass flux correction factor
    283367    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
    284368    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
     
    286370    REAL(wp) ::  time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information   
    287371    REAL(wp) ::  time_stg_call = 0.0_wp   !< time counter for calling generator   
     372   
     373    REAL(wp), DIMENSION(3) ::  mc_factor = 1.0_wp !< correction factor for the u,v,w-components to maintain original mass flux
    288374   
    289375   
     
    421507 SUBROUTINE stg_check_parameters
    422508
    423 
    424     USE control_parameters,                                                    &
    425         ONLY:  bc_lr, bc_ns, child_domain, nesting_offline, rans_mode,         &
    426                turbulent_inflow
    427 
    428     USE pmc_interface,                                                         &
    429         ONLY : rans_mode_parent
    430 
    431 
    432509    IMPLICIT NONE
    433510
     
    548625 SUBROUTINE stg_init
    549626
    550 
    551     USE arrays_3d,                                                             &
    552         ONLY:  ddzw, u_init, v_init
    553 
    554     USE control_parameters,                                                    &
    555         ONLY:  child_domain, coupling_char, e_init, nesting_offline, rans_mode
    556 
    557     USE grid_variables,                                                        &
    558         ONLY:  ddy
    559 
    560     USE indices,                                                               &
    561         ONLY:  nz
    562 
    563     USE pmc_interface,                                                         &
    564         ONLY : rans_mode_parent
    565 
    566 
    567627    IMPLICIT NONE
    568628
     
    581641    INTEGER(iwp) :: nseed                    !< dimension of random number seed
    582642    INTEGER(iwp) :: startseed = 1234567890   !< start seed for random number generator
    583 
     643   
     644    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_xz_l = 0 !< number of non-topography grid points at xz-cross-section on subdomain
     645    INTEGER(iwp), DIMENSION(3) ::  nr_non_topo_yz_l = 0 !< number of non-topography grid points at yz-cross-section on subdomain
    584646!
    585647!-- Dummy variables used for reading profiles
     
    782844    ELSE
    783845!
     846!--    Define length scale for the imposed turbulence, which is defined as
     847!--    8 times the minimum grid spacing
     848       length_scale = 8.0_wp * MIN( dx, dy, MINVAL( dzw ) )
     849!
     850!--    Define constant to gradually decrease length scales and Reynolds stress
     851!--    above ABL top. Assure that no zero length scales are used.
     852       d_l = blend_coeff / MAX( length_scale, dx, dy, MINVAL( dzw ) )
     853!
    784854!--    Set flag indicating that turbulence is parametrized
    785855       parametrize_inflow_turbulence = .TRUE.
     856!
     857!--    In case of dirichlet inflow boundary conditions only at one lateral
     858!--    boundary, i.e. in the case no offline or self nesting is applied but
     859!--    synthetic turbulence shall be parametrized nevertheless, the
     860!--    boundary-layer depth need to determined first.
     861       IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )                &
     862          CALL nesting_offl_calc_zi
    786863!
    787864!--    Determine boundary-layer depth, which is used to initialize lenght
     
    897974             DO  k = nzb, nzt+1
    898975                IF  ( a11(k) .NE. 0._wp ) THEN
    899                    fu_yz(k,j) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
     976                   fu_yz(k,j) = ( u(k,j,i) - u_init(k) ) / a11(k)
    900977                ELSE
    901978                   fu_yz(k,j) = 0._wp
     
    903980
    904981                IF  ( a22(k) .NE. 0._wp ) THEN
    905                    fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) - &
    906                                v_init(k) ) / a22(k)
     982                   fv_yz(k,j) = ( v(k,j,i) -                                  &
     983                                  a21(k) * fu_yz(k,j) - v_init(k) ) / a22(k)
    907984                ELSE
    908985                   fv_yz(k,j) = 0._wp
     
    910987
    911988                IF  ( a33(k) .NE. 0._wp ) THEN
    912                    fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - &
    913                                a32(k) * fv_yz(k,j) ) / a33(k)
     989                   fw_yz(k,j) = ( w(k,j,i) -                                   &
     990                                  a31(k) * fu_yz(k,j) - a32(k) *               &
     991                                  fv_yz(k,j) ) / a33(k)
    914992                ELSE
    915993                   fw_yz = 0._wp
     
    9291007
    9301008                IF  ( a11(k) .NE. 0._wp ) THEN
    931                    fu_xz(k,i) = ( u(k,j,i) / mc_factor - u_init(k) ) / a11(k)
     1009                   fu_xz(k,i) = ( u(k,j,i) - u_init(k) ) / a11(k)
    9321010                ELSE
    9331011                   fu_xz(k,i) = 0._wp
     
    9351013
    9361014                IF  ( a22(k) .NE. 0._wp ) THEN
    937                    fv_xz(k,i) = ( v(k,j,i) / mc_factor - a21(k) * fu_xz(k,i) - &
    938                                v_init(k) ) / a22(k)
     1015                   fv_xz(k,i) = ( v(k,j,i) -                                  &
     1016                                  a21(k) * fu_xz(k,i) - v_init(k) ) / a22(k)
    9391017                ELSE
    9401018                   fv_xz(k,i) = 0._wp
     
    9421020
    9431021                IF  ( a33(k) .NE. 0._wp ) THEN
    944                    fw_xz(k,i) = ( w(k,j,i) / mc_factor - a31(k) * fu_xz(k,i) - &
    945                                a32(k) * fv_xz(k,i) ) / a33(k)
     1022                   fw_xz(k,i) = ( w(k,j,i) -                                   &
     1023                                  a31(k) * fu_xz(k,i) -                        &
     1024                                  a32(k) * fv_xz(k,i) ) / a33(k)
    9461025                ELSE
    9471026                   fw_xz = 0._wp
     
    9521031       ENDIF
    9531032    ENDIF
     1033!
     1034!-- Count the number of non-topography grid points at the boundaries where
     1035!-- perturbations are imposed. This number is later used for bias corrections
     1036!-- of the perturbations, i.e. to force that their mean is zero. Please note,
     1037!-- due to the asymetry of u and v along x and y direction, respectively,
     1038!-- different cases must be distinguished.
     1039    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
     1040!
     1041!--    Number of grid points where perturbations are imposed on u
     1042       IF ( myidx == id_stg_left  )  i = nxl
     1043       IF ( myidx == id_stg_right )  i = nxr+1
     1044       
     1045       nr_non_topo_yz_l(1) = SUM(                                              &
     1046                          MERGE( 1, 0,                                         &
     1047                                 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 1 ) ) )
     1048!
     1049!--    Number of grid points where perturbations are imposed on v and w                                   
     1050       IF ( myidx == id_stg_left  )  i = nxl-1
     1051       IF ( myidx == id_stg_right )  i = nxr+1
     1052       
     1053       nr_non_topo_yz_l(2) = SUM(                                              &
     1054                          MERGE( 1, 0,                                         &
     1055                                 BTEST( wall_flags_0(nzb:nzt,nysv:nyn,i), 2 ) ) )
     1056       nr_non_topo_yz_l(3) = SUM(                                              &
     1057                          MERGE( 1, 0,                                         &
     1058                                 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 3 ) ) )
     1059                                   
     1060#if defined( __parallel )
     1061       CALL MPI_ALLREDUCE( nr_non_topo_yz_l, nr_non_topo_yz, 3, MPI_INTEGER,   &
     1062                           MPI_SUM, comm1dy, ierr )
     1063#else
     1064       nr_non_topo_yz = nr_non_topo_yz_l
     1065#endif 
     1066    ENDIF
     1067   
     1068    IF ( myidy == id_stg_south  .OR.  myidy == id_stg_north )  THEN
     1069!
     1070!--    Number of grid points where perturbations are imposed on v
     1071       IF ( myidy == id_stg_south )  j = nys
     1072       IF ( myidy == id_stg_north )  j = nyn+1
     1073       
     1074       nr_non_topo_xz_l(2) = SUM(                                              &
     1075                          MERGE( 1, 0,                                         &
     1076                                 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
     1077!
     1078!--    Number of grid points where perturbations are imposed on u and w
     1079       IF ( myidy == id_stg_south )  j = nys-1
     1080       IF ( myidy == id_stg_north )  j = nyn+1
     1081       
     1082       nr_non_topo_xz_l(1) = SUM(                                              &
     1083                          MERGE( 1, 0,                                         &
     1084                                 BTEST( wall_flags_0(nzb:nzt,j,nxlu:nxr), 1 ) ) )
     1085       nr_non_topo_xz_l(3) = SUM(                                              &
     1086                          MERGE( 1, 0,                                         &
     1087                                 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 3 ) ) )
     1088                                   
     1089#if defined( __parallel )
     1090       CALL MPI_ALLREDUCE( nr_non_topo_xz_l, nr_non_topo_xz, 3, MPI_INTEGER,   &
     1091                           MPI_SUM, comm1dx, ierr )
     1092#else
     1093       nr_non_topo_xz = nr_non_topo_xz_l
     1094#endif 
     1095    ENDIF
     1096
    9541097
    9551098 END SUBROUTINE stg_init
     
    10591202
    10601203
    1061     USE control_parameters,                                                    &
    1062         ONLY: length, restart_string
    1063 
    1064 
    10651204    IMPLICIT NONE
    10661205
     
    10721211
    10731212    SELECT CASE ( restart_string(1:length) )
    1074 
    1075        CASE ( 'mc_factor' )
    1076           READ ( 13 )  mc_factor
     1213         
     1214       CASE ( 'time_stg_adjust' )
     1215          READ ( 13 )  time_stg_adjust
     1216         
     1217       CASE ( 'time_stg_call' )
     1218          READ ( 13 )  time_stg_call
     1219         
    10771220       CASE ( 'use_syn_turb_gen' )
    10781221          READ ( 13 )  use_syn_turb_gen
     
    10971240
    10981241    IMPLICIT NONE
    1099 
    1100     CALL wrd_write_string( 'mc_factor' )
    1101     WRITE ( 14 )  mc_factor
     1242   
     1243    CALL wrd_write_string( 'time_stg_adjust' )
     1244    WRITE ( 14 )  time_stg_adjust
     1245   
     1246    CALL wrd_write_string( 'time_stg_call' )
     1247    WRITE ( 14 )  time_stg_call
    11021248
    11031249    CALL wrd_write_string( 'use_syn_turb_gen' )
     
    11181264 SUBROUTINE stg_main
    11191265
    1120 
    1121     USE arrays_3d,                                                             &
    1122         ONLY:  dzw
    1123 
    1124     USE control_parameters,                                                    &
    1125         ONLY:  child_domain, dt_3d,                                            &
    1126                nesting_offline, rans_mode, time_since_reference_point,         &
    1127                volume_flow_initial
    1128 
    1129     USE grid_variables,                                                        &
    1130         ONLY:  dx, dy
    1131 
    1132     USE indices,                                                               &
    1133         ONLY:  wall_flags_0
    1134 
    1135     USE pmc_interface,                                                         &
    1136         ONLY : rans_mode_parent
    1137 
    1138 
    11391266    IMPLICIT NONE
    11401267
     
    11421269    INTEGER(iwp) :: j           !< loop index in y-direction
    11431270    INTEGER(iwp) :: k           !< loop index in z-direction
     1271   
     1272    LOGICAL  :: stg_call = .FALSE. !< control flag indicating whether turbulence was updated or only restored from last call
    11441273
    11451274    REAL(wp) :: dt_stg          !< wheighted subtimestep
    1146     REAL(wp) :: mc_factor_l     !< local mass flux correction factor
    1147     REAL(wp) :: volume_flow     !< mass flux through lateral boundary
    1148     REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
    1149 
     1275   
     1276    REAL(wp), DIMENSION(3) :: mc_factor_l   !< local mass flux correction factor
    11501277
    11511278    IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'start' )
     
    11531280!-- Calculate time step which is needed for filter functions
    11541281    dt_stg = MAX( dt_3d, dt_stg_call )
     1282!
     1283!-- Check if synthetic turbulence generator needs to be called and new
     1284!-- perturbations need to be created or if old disturbances can be imposed
     1285!-- again.
     1286    IF ( time_stg_call >= dt_stg_call  .AND.                                  &
     1287         intermediate_timestep_count == intermediate_timestep_count_max  )  THEN
     1288       stg_call = .TRUE.
     1289    ELSE
     1290       stg_call = .FALSE.
     1291    ENDIF
    11551292!
    11561293!-- Initial value of fu, fv, fw
     
    11831320!
    11841321!-- New set of fu, fv, fw
    1185     CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
    1186     CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
    1187     CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
    1188    
    1189     IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
    1190                           .AND.  .NOT.  rans_mode ) )  THEN
    1191 !
    1192 !--    Generate turbulence at right boundary
    1193        CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
    1194        CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
    1195        CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
    1196 !
    1197 !--    Generate turbulence at north boundary
    1198        CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
    1199        CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
    1200        CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
    1201 !
    1202 !--    Generate turbulence at south boundary
    1203        CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
    1204        CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
    1205        CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     1322    IF ( stg_call )  THEN
     1323!     
     1324!--    Generate turbulence at left boundary (required in all applications,
     1325!--    i.e. offline nesting and turbulent inflow from the left)
     1326       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left )
     1327       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
     1328       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
     1329       
     1330       IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent      &
     1331                             .AND.  .NOT.  rans_mode ) )  THEN
     1332!     
     1333!--       Generate turbulence at right boundary
     1334          CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
     1335          CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
     1336          CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
     1337!     
     1338!--       Generate turbulence at north boundary
     1339          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
     1340          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
     1341          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
     1342!     
     1343!--       Generate turbulence at south boundary
     1344          CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
     1345          CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
     1346          CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     1347       ENDIF
    12061348    ENDIF
    12071349   
     
    12091351!-- Turbulence generation at left and or right boundary
    12101352    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
    1211    
    1212        DO  j = nysg, nyng
    1213           DO  k = nzb, nzt + 1
    1214 !
    1215 !--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
    1216              IF ( tu(k) == 0.0_wp )  THEN
    1217                 fu_yz(k,j) = fuo_yz(k,j)
    1218              ELSE
    1219                 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
    1220                          fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
    1221              ENDIF
    1222 
    1223              IF ( tv(k) == 0.0_wp )  THEN
    1224                 fv_yz(k,j) = fvo_yz(k,j)
    1225              ELSE
    1226                 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
    1227                          fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    1228              ENDIF
    1229 
    1230              IF ( tw(k) == 0.0_wp )  THEN
    1231                 fw_yz(k,j) = fwo_yz(k,j)
    1232              ELSE
    1233                 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
    1234                          fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    1235              ENDIF
    1236 !
    1237 !--          Lund rotation following Eq. 17 in Xie and Castro (2008).
    1238 !--          Additional factors are added to improve the variance of v and w
    1239              IF( k == 0 )  THEN
    1240                 dist_yz(k,j,1) = 0.0_wp
    1241                 dist_yz(k,j,2) = 0.0_wp
    1242                 dist_yz(k,j,3) = 0.0_wp
    1243              ELSE
    1244                 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
    1245                 !experimental test of 1.2
     1353!
     1354!--    Calculate new set of perturbations. Do this only at last RK3-substep and
     1355!--    when dt_stg_call is exceeded. Else the old set of perturbations is
     1356!--    imposed
     1357       IF ( stg_call )  THEN
     1358       
     1359          DO  j = nysg, nyng
     1360             DO  k = nzb, nzt + 1
     1361!     
     1362!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
     1363                IF ( tu(k) == 0.0_wp )  THEN
     1364                   fu_yz(k,j) = fuo_yz(k,j)
     1365                ELSE
     1366                   fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
     1367                            fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
     1368                ENDIF
     1369       
     1370                IF ( tv(k) == 0.0_wp )  THEN
     1371                   fv_yz(k,j) = fvo_yz(k,j)
     1372                ELSE
     1373                   fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
     1374                            fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
     1375                ENDIF
     1376       
     1377                IF ( tw(k) == 0.0_wp )  THEN
     1378                   fw_yz(k,j) = fwo_yz(k,j)
     1379                ELSE
     1380                   fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
     1381                            fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
     1382                ENDIF
     1383             ENDDO
     1384          ENDDO
     1385         
     1386          dist_yz(nzb,:,1) = 0.0_wp
     1387          dist_yz(nzb,:,2) = 0.0_wp
     1388          dist_yz(nzb,:,3) = 0.0_wp
     1389         
     1390          IF ( myidx == id_stg_left  )  i = nxl
     1391          IF ( myidx == id_stg_right )  i = nxr+1       
     1392          DO  j = nysg, nyng
     1393             DO  k = nzb+1, nzt + 1
     1394!     
     1395!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
     1396!--             Additional factors are added to improve the variance of v and w
     1397                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) *          &
     1398                                    MERGE( 1.0_wp, 0.0_wp,                     &
     1399                                           BTEST( wall_flags_0(k,j,i), 1 ) )   
     1400             ENDDO
     1401          ENDDO
     1402         
     1403          IF ( myidx == id_stg_left  )  i = nxl-1
     1404          IF ( myidx == id_stg_right )  i = nxr+1
     1405          DO  j = nysg, nyng
     1406             DO  k = nzb+1, nzt + 1
     1407!     
     1408!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
     1409!--             Additional factors are added to improve the variance of v and w
     1410!--             experimental test of 1.2                                       
    12461411                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
    12471412                                      * 1.2_wp )                               &
    12481413                                     * (   a21(k) * fu_yz(k,j)                 &
    1249                                          + a22(k) * fv_yz(k,j) ), 3.0_wp )
     1414                                         + a22(k) * fv_yz(k,j) ), 3.0_wp ) *   &
     1415                                    MERGE( 1.0_wp, 0.0_wp,                     &
     1416                                           BTEST( wall_flags_0(k,j,i), 2 ) )   
    12501417                dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
    12511418                                      * 1.3_wp )                               &
    12521419                                    * (   a31(k) * fu_yz(k,j)                  &
    12531420                                        + a32(k) * fv_yz(k,j)                  &
    1254                                         + a33(k) * fw_yz(k,j) ), 3.0_wp )
    1255              ENDIF
    1256 
     1421                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )  *   &
     1422                                    MERGE( 1.0_wp, 0.0_wp,                     &
     1423                                           BTEST( wall_flags_0(k,j,i), 3 ) )
     1424             ENDDO
    12571425          ENDDO
    1258        ENDDO
     1426       ENDIF
    12591427
    12601428!
    12611429!--    Mass flux correction following Kim et al. (2013)
    12621430!--    This correction factor insures that the mass flux is preserved at the
    1263 !--    inflow boundary
    1264        IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN
    1265           mc_factor_l = 0.0_wp
    1266           mc_factor   = 0.0_wp
    1267           DO  j = nys, nyn
    1268              DO  k = nzb+1, nzt
    1269                 mc_factor_l = mc_factor_l + dzw(k)  *                          &
    1270                               ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) )
    1271              ENDDO
    1272           ENDDO
    1273 
    1274 #if defined( __parallel )
    1275           CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
    1276                               comm1dy, ierr )
    1277 #else
    1278           mc_factor = mc_factor_l
    1279 #endif
    1280 
    1281           mc_factor = volume_flow_initial(1) / mc_factor
    1282 
    1283 !
    1284 !--       Add disturbance at the inflow
    1285           DO  j = nysg, nyng
    1286              DO  k = nzb, nzt+1
    1287                  u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +              &
    1288                                       dist_yz(k,j,1)             ) * mc_factor
    1289                  v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +              &
    1290                                       dist_yz(k,j,2)             ) * mc_factor
    1291                  w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)               * mc_factor
    1292              ENDDO
    1293           ENDDO
    1294 
    1295        ELSE
    1296 !
    1297 !--       First, calculate volume flow at yz boundary
     1431!--    inflow boundary. First, calculate mean value of the imposed
     1432!--    perturbations at yz boundary.
     1433!--    Note, this needs to be done only after the last RK3-substep, else
     1434!--    the boundary values won't be accessed.
     1435       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
     1436          mc_factor_l   = 0.0_wp
     1437          mc_factor     = 0.0_wp
     1438!       
     1439!--       Sum up the original volume flows (with and without perturbations).
     1440!--       Note, for non-normal components (here v and w) it is actually no
     1441!--       volume flow.
    12981442          IF ( myidx == id_stg_left  )  i = nxl
    12991443          IF ( myidx == id_stg_right )  i = nxr+1
    1300 
    1301           volume_flow_l = 0.0_wp
    1302           volume_flow   = 0.0_wp
    1303           mc_factor_l   = 0.0_wp
    1304           mc_factor     = 0.0_wp
    1305           DO  j = nys, nyn
    1306              DO  k = nzb+1, nzt
    1307                 volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy         &
    1308                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    1309                                               BTEST( wall_flags_0(k,j,i), 1 ) )
    1310 
    1311                 mc_factor_l = mc_factor_l     + ( u(k,j,i) + dist_yz(k,j,1) )  &
    1312                                                          * dzw(k) * dy         &
    1313                                      * MERGE( 1.0_wp, 0.0_wp,                  &
    1314                                               BTEST( wall_flags_0(k,j,i), 1 ) )
    1315              ENDDO
    1316           ENDDO
     1444         
     1445          mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1)                     &
     1446                         * MERGE( 1.0_wp, 0.0_wp,                              &
     1447                                  BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 1 ) ) )
     1448       
     1449          IF ( myidx == id_stg_left  )  i = nxl-1
     1450          IF ( myidx == id_stg_right )  i = nxr+1
     1451         
     1452          mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2)                    &
     1453                          * MERGE( 1.0_wp, 0.0_wp,                             &
     1454                                   BTEST( wall_flags_0(nzb:nzt,nysv:nyn,i), 2 ) ) )
     1455          mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3)                     &
     1456                          * MERGE( 1.0_wp, 0.0_wp,                             &
     1457                                   BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 3 ) ) )
     1458         
    13171459#if defined( __parallel )
    1318           CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                      &
    1319                               1, MPI_REAL, MPI_SUM, comm1dy, ierr )
    13201460          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
    1321                               1, MPI_REAL, MPI_SUM, comm1dy, ierr )
    1322 #else
    1323           volume_flow = volume_flow_l
    1324           mc_factor   = mc_factor_l
     1461                              3, MPI_REAL, MPI_SUM, comm1dy, ierr )           
     1462#else                                                                         
     1463          mc_factor   = mc_factor_l                                           
    13251464#endif
    1326 
    1327           mc_factor = volume_flow / mc_factor
    1328 
    1329 !
     1465!     
     1466!--       Calculate correction factor and force zero mean perturbations.
     1467          mc_factor = mc_factor / REAL( nr_non_topo_yz, KIND = wp )           
     1468                                                                               
     1469          IF ( myidx == id_stg_left  )  i = nxl                               
     1470          IF ( myidx == id_stg_right )  i = nxr+1                             
     1471                                                                               
     1472          dist_yz(:,:,1) = ( dist_yz(:,:,1) - mc_factor(1) )                   &
     1473                        * MERGE( 1.0_wp, 0.0_wp,                               &
     1474                                 BTEST( wall_flags_0(:,:,i), 1 ) )             
     1475                                                                               
     1476                                                                               
     1477          IF ( myidx == id_stg_left  )  i = nxl-1                             
     1478          IF ( myidx == id_stg_right )  i = nxr+1                             
     1479                                                                               
     1480          dist_yz(:,:,2) = ( dist_yz(:,:,2) - mc_factor(2) )                   &
     1481                        * MERGE( 1.0_wp, 0.0_wp,                               &
     1482                                 BTEST( wall_flags_0(:,:,i), 2 ) )             
     1483                                                                               
     1484          dist_yz(:,:,3) = ( dist_yz(:,:,3) - mc_factor(3) )                   &
     1485                        * MERGE( 1.0_wp, 0.0_wp,                               &
     1486                                 BTEST( wall_flags_0(:,:,i), 3 ) )
     1487!     
    13301488!--       Add disturbances
    13311489          IF ( myidx == id_stg_left  )  THEN
    1332              DO  j = nys, nyn
    1333                 DO  k = nzb+1, nzt
    1334                    u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
    1335                        * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
    1336                                             BTEST( wall_flags_0(k,j,0), 1 ) )
    1337                    u(k,j,-1) = u(k,j,0)
    1338                    v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
    1339                        * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
    1340                                             BTEST( wall_flags_0(k,j,-1), 2 ) )
    1341                    w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3)  )               &
    1342                        * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
     1490!     
     1491!--          For the left boundary distinguish between mesoscale offline / self
     1492!--          nesting and turbulent inflow at the left boundary only. In the latter
     1493!--          case turbulence is imposed on the mean inflow profiles.
     1494             IF ( .NOT. nesting_offline  .AND.  .NOT. child_domain )  THEN
     1495!           
     1496!--             Add disturbance at the inflow
     1497                DO  j = nysg, nyng
     1498                   DO  k = nzb, nzt+1
     1499                      u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) +         &
     1500                                           dist_yz(k,j,1)             )        &
     1501                                   * MERGE( 1.0_wp, 0.0_wp,                    &
     1502                                            BTEST( wall_flags_0(k,j,0), 1 ) ) 
     1503                      v(k,j,-nbgp:-1)  = ( mean_inflow_profiles(k,2) +         &
     1504                                           dist_yz(k,j,2)             )        &
     1505                                   * MERGE( 1.0_wp, 0.0_wp,                    &
     1506                                            BTEST( wall_flags_0(k,j,-1), 2 ) )
     1507                      w(k,j,-nbgp:-1)  =   dist_yz(k,j,3)                      &
     1508                                   * MERGE( 1.0_wp, 0.0_wp,                    &
    13431509                                            BTEST( wall_flags_0(k,j,-1), 3 ) )
     1510                   ENDDO
    13441511                ENDDO
    1345              ENDDO
     1512             ELSE
     1513             
     1514                DO  j = nys, nyn
     1515                   DO  k = nzb+1, nzt
     1516                      u(k,j,0)   = ( u(k,j,0) + dist_yz(k,j,1) )               &
     1517                                 * MERGE( 1.0_wp, 0.0_wp,                      &
     1518                                          BTEST( wall_flags_0(k,j,0), 1 ) )
     1519                      u(k,j,-1)  = u(k,j,0)
     1520                      v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )            &
     1521                                 * MERGE( 1.0_wp, 0.0_wp,                      &
     1522                                          BTEST( wall_flags_0(k,j,-1), 2 ) )
     1523                      w(k,j,-1)  = ( w(k,j,-1)  + dist_yz(k,j,3) )             &
     1524                                 * MERGE( 1.0_wp, 0.0_wp,                      &
     1525                                          BTEST( wall_flags_0(k,j,-1), 3 ) )
     1526                   ENDDO
     1527                ENDDO
     1528             ENDIF
    13461529          ENDIF
    13471530          IF ( myidx == id_stg_right  )  THEN
     
    13491532                DO  k = nzb+1, nzt
    13501533                   u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
    1351                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    1352                                           BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
     1534                                  * MERGE( 1.0_wp, 0.0_wp,                     &
     1535                                           BTEST( wall_flags_0(k,j,nxr+1), 1 ) )
    13531536                   v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) )            &
    1354                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    1355                                           BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
     1537                                  * MERGE( 1.0_wp, 0.0_wp,                     &
     1538                                           BTEST( wall_flags_0(k,j,nxr+1), 2 ) )
    13561539                   w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) )            &
    1357                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    1358                                           BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
     1540                                  * MERGE( 1.0_wp, 0.0_wp,                     &
     1541                                           BTEST( wall_flags_0(k,j,nxr+1), 3 ) )
    13591542                ENDDO
    13601543             ENDDO
    13611544          ENDIF
    13621545       ENDIF
    1363 
    13641546    ENDIF
    13651547!
    13661548!-- Turbulence generation at north and south boundary
    13671549    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
    1368    
    1369        DO  i = nxlg, nxrg
    1370           DO  k = nzb, nzt + 1
    1371 !
    1372 !--          Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
    1373              IF ( tu(k) == 0.0_wp )  THEN
    1374                 fu_xz(k,i) = fuo_xz(k,i)
    1375              ELSE
    1376                 fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
    1377                          fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
    1378              ENDIF
    1379 
    1380              IF ( tv(k) == 0.0_wp )  THEN
    1381                 fv_xz(k,i) = fvo_xz(k,i)
    1382              ELSE
    1383                 fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
    1384                       fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    1385              ENDIF
    1386            
    1387              IF ( tw(k) == 0.0_wp )  THEN
    1388                 fw_xz(k,i) = fwo_xz(k,i)
    1389              ELSE
    1390                 fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
    1391                       fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    1392              ENDIF
    1393 !
    1394 !--          Lund rotation following Eq. 17 in Xie and Castro (2008).
    1395 !--          Additional factors are added to improve the variance of v and w
    1396              IF( k == 0 )  THEN
    1397                 dist_xz(k,i,1) = 0.0_wp
    1398                 dist_xz(k,i,2) = 0.0_wp
    1399                 dist_xz(k,i,3) = 0.0_wp
    1400 
    1401              ELSE
    1402                 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
    1403                 !experimental test of 1.2
     1550!
     1551!--    Calculate new set of perturbations. Do this only at last RK3-substep and
     1552!--    when dt_stg_call is exceeded. Else the old set of perturbations is
     1553!--    imposed
     1554       IF ( stg_call )  THEN
     1555          DO  i = nxlg, nxrg
     1556             DO  k = nzb, nzt + 1
     1557!         
     1558!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
     1559                IF ( tu(k) == 0.0_wp )  THEN
     1560                   fu_xz(k,i) = fuo_xz(k,i)
     1561                ELSE
     1562                   fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
     1563                            fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
     1564                ENDIF
     1565         
     1566                IF ( tv(k) == 0.0_wp )  THEN
     1567                   fv_xz(k,i) = fvo_xz(k,i)
     1568                ELSE
     1569                   fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
     1570                         fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
     1571                ENDIF
     1572               
     1573                IF ( tw(k) == 0.0_wp )  THEN
     1574                   fw_xz(k,i) = fwo_xz(k,i)
     1575                ELSE
     1576                   fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
     1577                         fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
     1578                ENDIF
     1579             ENDDO
     1580          ENDDO
     1581         
     1582         
     1583          dist_xz(nzb,:,1) = 0.0_wp
     1584          dist_xz(nzb,:,2) = 0.0_wp
     1585          dist_xz(nzb,:,3) = 0.0_wp
     1586         
     1587          IF ( myidy == id_stg_south  ) j = nys
     1588          IF ( myidy == id_stg_north )  j = nyn+1
     1589          DO  i = nxlg, nxrg
     1590             DO  k = nzb+1, nzt + 1
     1591!         
     1592!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
     1593!--             Additional factors are added to improve the variance of v and w
     1594                !experimental test of 1.2                                         
    14041595                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
    14051596                                      * 1.2_wp )                               &
    14061597                                     * (   a21(k) * fu_xz(k,i)                 &
    1407                                          + a22(k) * fv_xz(k,i) ), 3.0_wp )
     1598                                         + a22(k) * fv_xz(k,i) ), 3.0_wp ) *   &
     1599                                    MERGE( 1.0_wp, 0.0_wp,                     &
     1600                                           BTEST( wall_flags_0(k,j,i), 2 ) )   
     1601             ENDDO
     1602          ENDDO
     1603         
     1604          IF ( myidy == id_stg_south  ) j = nys-1
     1605          IF ( myidy == id_stg_north )  j = nyn+1
     1606          DO  i = nxlg, nxrg
     1607             DO  k = nzb+1, nzt + 1
     1608!         
     1609!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
     1610!--             Additional factors are added to improve the variance of v and w
     1611                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) *          &
     1612                                    MERGE( 1.0_wp, 0.0_wp,                     &
     1613                                           BTEST( wall_flags_0(k,j,i), 1 ) )   
     1614         
    14081615                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
    14091616                                      * 1.3_wp )                               &
    14101617                                    * (   a31(k) * fu_xz(k,i)                  &
    14111618                                        + a32(k) * fv_xz(k,i)                  &
    1412                                         + a33(k) * fw_xz(k,i) ), 3.0_wp )
    1413              ENDIF
    1414 
    1415           ENDDO
    1416        ENDDO
    1417 
    1418 !
    1419 !--    Mass flux correction following Kim et al. (2013)
    1420 !--    This correction factor insures that the mass flux is preserved at the
    1421 !--    inflow boundary.
    1422 !--    First, calculate volume flow at xz boundary
    1423        IF ( myidy == id_stg_south  ) j = nys
    1424        IF ( myidy == id_stg_north )  j = nyn+1
    1425 
    1426        volume_flow_l = 0.0_wp
    1427        volume_flow   = 0.0_wp
    1428        mc_factor_l   = 0.0_wp
    1429        mc_factor     = 0.0_wp
    1430        DO  i = nxl, nxr
    1431           DO  k = nzb+1, nzt
    1432              volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx            &
    1433                                   * MERGE( 1.0_wp, 0.0_wp,                     &
    1434                                            BTEST( wall_flags_0(k,j,i), 2 ) )
    1435 
    1436              mc_factor_l = mc_factor_l     + ( v(k,j,i) + dist_xz(k,i,2) )     &
    1437                                                       * dzw(k) * dx            &
    1438                                   * MERGE( 1.0_wp, 0.0_wp,                     &
    1439                                            BTEST( wall_flags_0(k,j,i), 2 ) )
    1440           ENDDO
    1441        ENDDO
    1442 #if defined( __parallel )
    1443        CALL MPI_ALLREDUCE( volume_flow_l, volume_flow,                         &
    1444                            1, MPI_REAL, MPI_SUM, comm1dx, ierr )
    1445        CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                             &
    1446                            1, MPI_REAL, MPI_SUM, comm1dx, ierr )
    1447 #else
    1448        volume_flow = volume_flow_l
    1449        mc_factor   = mc_factor_l
    1450 #endif
    1451 
    1452        mc_factor = volume_flow / mc_factor
    1453 
    1454 !
    1455 !--    Add disturbances
    1456        IF ( myidy == id_stg_south  )  THEN
    1457 
    1458           DO  i = nxl, nxr
    1459              DO  k = nzb+1, nzt
    1460                 u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
    1461                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
    1462                                            BTEST( wall_flags_0(k,-1,i), 1 ) )
    1463                 v(k,0,i)  = ( v(k,0,i)  + dist_xz(k,i,2)  )                    &
    1464                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
    1465                                            BTEST( wall_flags_0(k,0,i), 2 ) )
    1466                 v(k,-1,i) = v(k,0,i)
    1467                 w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
    1468                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
    1469                                            BTEST( wall_flags_0(k,-1,i), 3 ) )
     1619                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )  *   &
     1620                                    MERGE( 1.0_wp, 0.0_wp,                     &
     1621                                           BTEST( wall_flags_0(k,j,i), 3 ) )
    14701622             ENDDO
    14711623          ENDDO
    14721624       ENDIF
    1473        IF ( myidy == id_stg_north  )  THEN
    1474 
    1475           DO  i = nxl, nxr
    1476              DO  k = nzb+1, nzt
    1477                 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
    1478                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    1479                                           BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
    1480                 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
    1481                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    1482                                           BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
    1483                 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
    1484                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    1485                                           BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
     1625
     1626!
     1627!--    Mass flux correction following Kim et al. (2013)
     1628!--    This correction factor insures that the mass flux is preserved at the
     1629!--    inflow boundary. First, calculate mean value of the imposed
     1630!--    perturbations at yz boundary.
     1631!--    Note, this needs to be done only after the last RK3-substep, else
     1632!--    the boundary values won't be accessed.
     1633       IF ( intermediate_timestep_count == intermediate_timestep_count_max )  THEN
     1634          mc_factor_l   = 0.0_wp
     1635          mc_factor     = 0.0_wp
     1636         
     1637          IF ( myidy == id_stg_south  ) j = nys
     1638          IF ( myidy == id_stg_north )  j = nyn+1
     1639         
     1640          mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2)                     &
     1641                         * MERGE( 1.0_wp, 0.0_wp,                              &
     1642                                  BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 2 ) ) )
     1643         
     1644          IF ( myidy == id_stg_south  ) j = nys-1
     1645          IF ( myidy == id_stg_north )  j = nyn+1
     1646         
     1647          mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1)                    &
     1648                         * MERGE( 1.0_wp, 0.0_wp,                              &
     1649                                  BTEST( wall_flags_0(nzb:nzt,j,nxlu:nxr), 1 ) ) )
     1650          mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3)                     &
     1651                         * MERGE( 1.0_wp, 0.0_wp,                              &
     1652                                  BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 3 ) ) )
     1653         
     1654#if defined( __parallel )
     1655          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor,                          &
     1656                              3, MPI_REAL, MPI_SUM, comm1dx, ierr )
     1657#else     
     1658          mc_factor   = mc_factor_l
     1659#endif
     1660         
     1661          mc_factor = mc_factor / REAL( nr_non_topo_xz, KIND = wp )
     1662         
     1663          IF ( myidy == id_stg_south  ) j = nys
     1664          IF ( myidy == id_stg_north )  j = nyn+1
     1665         
     1666          dist_xz(:,:,2)   = ( dist_xz(:,:,2) - mc_factor(2) )                 &
     1667                           * MERGE( 1.0_wp, 0.0_wp,                            &
     1668                                    BTEST( wall_flags_0(:,j,:), 2 ) )         
     1669                                                                               
     1670                                                                               
     1671          IF ( myidy == id_stg_south  ) j = nys-1                             
     1672          IF ( myidy == id_stg_north )  j = nyn+1                             
     1673                                                                               
     1674          dist_xz(:,:,1)   = ( dist_xz(:,:,1) - mc_factor(1) )                 &
     1675                           * MERGE( 1.0_wp, 0.0_wp,                            &
     1676                                    BTEST( wall_flags_0(:,j,:), 1 ) )         
     1677                                                                               
     1678          dist_xz(:,:,3)   = ( dist_xz(:,:,3) - mc_factor(3) )                 &
     1679                           * MERGE( 1.0_wp, 0.0_wp,                            &
     1680                                    BTEST( wall_flags_0(:,j,:), 3 ) )     
     1681!         
     1682!--       Add disturbances
     1683          IF ( myidy == id_stg_south  )  THEN
     1684             DO  i = nxl, nxr
     1685                DO  k = nzb+1, nzt
     1686                   u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                  &
     1687                                        * MERGE( 1.0_wp, 0.0_wp,               &
     1688                                                 BTEST( wall_flags_0(k,-1,i), 1 ) )
     1689                   v(k,0,i)  = ( v(k,0,i)  + dist_xz(k,i,2)  )                 &
     1690                                        * MERGE( 1.0_wp, 0.0_wp,               &
     1691                                                 BTEST( wall_flags_0(k,0,i), 2 ) )
     1692                   v(k,-1,i) = v(k,0,i)
     1693                   w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                 &
     1694                                        * MERGE( 1.0_wp, 0.0_wp,               &
     1695                                                 BTEST( wall_flags_0(k,-1,i), 3 ) )
     1696                ENDDO
    14861697             ENDDO
    1487           ENDDO
     1698          ENDIF
     1699          IF ( myidy == id_stg_north  )  THEN
     1700         
     1701             DO  i = nxl, nxr
     1702                DO  k = nzb+1, nzt
     1703                   u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )            &
     1704                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     1705                                             BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
     1706                   v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )            &
     1707                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     1708                                              BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
     1709                   w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )            &
     1710                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     1711                                              BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
     1712                ENDDO
     1713             ENDDO
     1714          ENDIF
    14881715       ENDIF
    14891716    ENDIF
    14901717!
    14911718!-- Finally, set time counter for calling STG to zero
    1492     time_stg_call = 0.0_wp
     1719    IF ( stg_call )  time_stg_call = 0.0_wp
    14931720
    14941721    IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'end' )
     
    15061733!------------------------------------------------------------------------------!
    15071734 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id )
    1508 
    1509 
    1510     USE indices,                                                               &
    1511         ONLY: ny
    15121735
    15131736    IMPLICIT NONE
     
    16331856 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id )
    16341857
    1635 
    1636     USE indices,                                                               &
    1637         ONLY: nx
    1638 
    1639 
    16401858    IMPLICIT NONE
    16411859
     
    17581976 SUBROUTINE parametrize_reynolds_stress
    17591977
    1760     USE arrays_3d,                                                             &
    1761         ONLY:  zu
    1762 
    17631978    IMPLICIT NONE
    17641979
     
    17671982    REAL(wp)     ::  zzi         !< ratio of z/zi
    17681983   
    1769 !
    1770 !--
    17711984    DO  k = nzb+1, nzt+1
    1772          
    1773        IF ( zu(k) <= zi_ribulk )  THEN
    1774 !
    1775 !--       Determine normalized height coordinate
    1776           zzi = zu(k) / zi_ribulk
    1777 !
    1778 !--       u'u' and v'v'. Assume isotropy. Note, add a small negative number
    1779 !--       to the denominator, else the merge-function can crash if scale_l is
    1780 !--       zero.
    1781           r11(k) = scale_us**2 * (                                             &
    1782                       MERGE( 0.35_wp * ABS(                                    &
    1783                            - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
    1784                                           )**( 2.0_wp / 3.0_wp ),              &
    1785                              0.0_wp,                                           &
    1786                              scale_l < 0.0_wp )                                &
    1787                     + 5.0_wp - 4.0_wp * zzi                                    &
    1788                                  )
    1789                                  
    1790           r22(k) = r11(k)
    1791 !
    1792 !--       w'w'
    1793           r33(k) = scale_wm**2 * (                                             &
    1794                       1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) &
    1795                     + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2            &                     
    1796                                  )
    1797 !
    1798 !--       u'w' and v'w'. Assume isotropy.
    1799           r31(k) = - scale_us**2 * (                                           &
    1800                          1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
    1801                                    )
    1802                              
    1803           r32(k) = r31(k)
    1804 !
    1805 !--       For u'v' no parametrization exist so far - ?. For simplicity assume
    1806 !--       a similar profile as for u'w'.
    1807           r21(k) = r31(k)
    1808 !
    1809 !--    Above the boundary layer, assmume laminar flow conditions.
    1810        ELSE
    1811           r11(k) = 10E-8_wp
    1812           r22(k) = 10E-8_wp
    1813           r33(k) = 10E-8_wp
    1814           r21(k) = 10E-8_wp
    1815           r31(k) = 10E-8_wp
    1816           r32(k) = 10E-8_wp
    1817        ENDIF
     1985!
     1986!--    Calculate function to gradually decrease Reynolds stress above ABL top.
     1987!--    The function decreases to 1/10 after one length scale above the
     1988!--    ABL top.
     1989       blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) )
     1990!
     1991!--    Determine normalized height coordinate
     1992       zzi = zu(k) / zi_ribulk
     1993!
     1994!--    u'u' and v'v'. Assume isotropy. Note, add a small negative number
     1995!--    to the denominator, else the merge-function can crash if scale_l is
     1996!--    zero.
     1997       r11(k) = scale_us**2 * (                                                &
     1998                   MERGE( 0.35_wp * ABS(                                       &
     1999                        - zi_ribulk / ( kappa * scale_l - 10E-4_wp )           &
     2000                                       )**( 2.0_wp / 3.0_wp ),                 &
     2001                          0.0_wp,                                              &
     2002                          scale_l < 0.0_wp )                                   &
     2003                 + 5.0_wp - 4.0_wp * zzi                                       &
     2004                              ) * blend                                       
     2005                                                                               
     2006       r22(k) = r11(k)                                                         
     2007!                                                                             
     2008!--    w'w'                                                                   
     2009       r33(k) = scale_wm**2 * (                                                &
     2010                   1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi )    &
     2011                 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2               &                     
     2012                              )  * blend                                       
     2013!                                                                             
     2014!--    u'w' and v'w'. Assume isotropy.                                         
     2015       r31(k) = - scale_us**2 * ( 1.01_wp - MIN( zzi, 1.0_wp ) ) * blend
     2016       r32(k) = r31(k)
     2017!
     2018!--    For u'v' no parametrization exist so far - ?. For simplicity assume
     2019!--    a similar profile as for u'w'.
     2020       r21(k) = r31(k)
    18182021    ENDDO
    18192022
     
    18452048!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
    18462049    DO  k = nzb+1, nzt+1
    1847        IF ( r11(k) > 0.0_wp )  THEN
     2050       IF ( r11(k) > 10E-6_wp )  THEN
    18482051          a11(k) = SQRT( r11(k) )
    18492052          a21(k) = r21(k) / a11(k)
     
    18572060    DO  k = nzb+1, nzt+1
    18582061       a22(k) = r22(k) - a21(k)**2
    1859        IF ( a22(k) > 0.0_wp )  THEN
     2062       IF ( a22(k) > 10E-6_wp )  THEN
    18602063          a22(k) = SQRT( a22(k) )
    18612064          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
     
    18672070    DO  k = nzb+1, nzt+1
    18682071       a33(k) = r33(k) - a31(k)**2 - a32(k)**2
    1869        IF ( a33(k) > 0.0_wp )  THEN
     2072       IF ( a33(k) > 10E-6_wp )  THEN
    18702073          a33(k) =  SQRT( a33(k) )
    18712074       ELSE
     
    18972100    IF ( debug_output_timestep )  CALL debug_message( 'stg_adjust', 'start' )
    18982101!
    1899 !-- Compute mean boundary layer height according to Richardson-Bulk
    1900 !-- criterion using the inflow profiles. Further velocity scale as well as
    1901 !-- mean friction velocity are calculated.
     2102!-- In case of dirichlet inflow boundary conditions only at one lateral
     2103!-- boundary, i.e. in the case no offline or self nesting is applied but
     2104!-- synthetic turbulence shall be parametrized nevertheless, the
     2105!-- boundary-layer depth need to determined first.
     2106    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )                   &
     2107       CALL nesting_offl_calc_zi   
     2108!
     2109!-- Compute scaling parameters (domain-averaged), such as friction velocity
     2110!-- are calculated.
    19022111    CALL calc_scaling_variables
    19032112!
     
    19422151 SUBROUTINE calc_length_and_time_scale
    19432152
    1944     USE arrays_3d,                                                             &
    1945         ONLY:  dzw, ddzw, u_init, v_init
    1946 
    1947     USE grid_variables,                                                        &
    1948         ONLY:  ddx, ddy, dx, dy
    1949 
    19502153    IMPLICIT NONE
    19512154
    19522155
    1953     INTEGER(iwp) :: k            !< loop index in z-direction
    1954 
    1955     REAL(wp) ::  length_scale    !< typical length scale
     2156    INTEGER(iwp) ::  k !< loop index in z-direction
     2157   
     2158    REAL(wp) ::  length_scale_dum     !< effectively used length scale
    19562159   
    19572160!
     
    19702173!-- using the length scales and the mean profiles of u- and v-component.
    19712174    DO  k = nzb+1, nzt+1
    1972    
    1973        length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
    1974          
    1975        IF ( zu(k) <= zi_ribulk )  THEN 
    1976 !
    1977 !--       Assume isotropic turbulence length scales
    1978           nux(k) = MAX( INT( length_scale * ddx     ), 1 )
    1979           nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
    1980           nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
    1981           nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
    1982           nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
    1983           nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
    1984           nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
    1985           nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
    1986           nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
    1987 !
    1988 !--       Limit time scales, else they become very larger for low wind speed,
    1989 !--       imposing long-living inflow perturbations which in turn propagate
    1990 !--       further into the model domain. Use u_init and v_init to calculate
    1991 !--       the time scales, which will be equal to the inflow profiles, both,
    1992 !--       in offline nesting mode or in dirichlet/radiation mode.
    1993           tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
    1994                                   ( ABS( u_init(k) ) + 0.1_wp ) )
    1995           tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
    1996                                   ( ABS( v_init(k) ) + 0.1_wp ) )
    1997 !
    1998 !--       Time scale of w-component is a mixture from u- and v-component.
    1999           tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
    2000 !
    2001 !--    Above the boundary layer length scales are zero, i.e. imposed turbulence
    2002 !--    is not correlated in space and time, just white noise. This saves
    2003 !--    computations power.
    2004        ELSE
    2005           nux(k) = 0.0_wp
    2006           nuy(k) = 0.0_wp
    2007           nuz(k) = 0.0_wp
    2008           nvx(k) = 0.0_wp
    2009           nvy(k) = 0.0_wp
    2010           nvz(k) = 0.0_wp
    2011           nwx(k) = 0.0_wp
    2012           nwy(k) = 0.0_wp
    2013           nwz(k) = 0.0_wp
    2014          
    2015           tu(k)  = 0.0_wp
    2016           tv(k)  = 0.0_wp
    2017           tw(k)  = 0.0_wp
    2018        ENDIF
     2175!
     2176!--    Determine blending parameter. Within the boundary layer length scales
     2177!--    are constant, while above lengths scales approach gradully zero,
     2178!--    i.e. imposed turbulence is not correlated in space and time,
     2179!--    just white noise, which saves computations power as the loops for the
     2180!--    computation of the filter functions depend on the length scales.
     2181!--    The value decreases to 1/10 after one length scale above the
     2182!--    ABL top.
     2183       blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) )
     2184!
     2185!--    Assume isotropic turbulence length scales
     2186       nux(k) = MAX( INT( length_scale * ddx     ), 1 ) * blend
     2187       nuy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend       
     2188       nvx(k) = MAX( INT( length_scale * ddx     ), 1 ) * blend
     2189       nvy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend       
     2190       nwx(k) = MAX( INT( length_scale * ddx     ), 1 ) * blend
     2191       nwy(k) = MAX( INT( length_scale * ddy     ), 1 ) * blend       
     2192!
     2193!--    Along the vertical direction limit the length scale further by the
     2194!--    boundary-layer depth to assure that no length scales larger than
     2195!--    the boundary-layer depth are used
     2196       length_scale_dum = MIN( length_scale, zi_ribulk )
     2197       
     2198       nuz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
     2199       nvz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
     2200       nwz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
     2201!
     2202!--    Limit time scales, else they become very larger for low wind speed,
     2203!--    imposing long-living inflow perturbations which in turn propagate
     2204!--    further into the model domain. Use u_init and v_init to calculate
     2205!--    the time scales, which will be equal to the inflow profiles, both,
     2206!--    in offline nesting mode or in dirichlet/radiation mode.
     2207       tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
     2208                               ( ABS( u_init(k) ) + 0.1_wp ) ) * blend
     2209       tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
     2210                               ( ABS( v_init(k) ) + 0.1_wp ) ) * blend
     2211!
     2212!--    Time scale of w-component is a mixture from u- and v-component.
     2213       tw(k)  = SQRT( tu(k)**2 + tv(k)**2 ) * blend
     2214     
    20192215    ENDDO
    20202216!
     
    20452241!------------------------------------------------------------------------------!
    20462242 SUBROUTINE calc_scaling_variables
    2047              
    2048     USE arrays_3d,                                                          &
    2049         ONLY:  drho_air         
    2050              
    2051     USE indices,                                                               &
    2052         ONLY:  nx, ny 
    20532243       
    20542244    USE surface_mod,                                                           &
  • palm/trunk/SOURCE/time_integration.f90

    r4017 r4022  
    2525! -----------------
    2626! $Id$
     27! Call synthetic turbulence generator at last RK3 substep right after boundary
     28! conditions are updated in offline nesting in order to assure that
     29! perturbations are always imposed
     30!
     31! 4017 2019-06-06 12:16:46Z schwenkel
    2732! Mass (volume) flux correction included to ensure global mass conservation for child domains.
    2833!
     
    12401245             CALL nesting_offl_bc
    12411246!
    1242 !--       Impose a turbulent inflow using synthetic generated turbulence,
    1243 !--       only once per time step.
    1244           IF ( use_syn_turb_gen  .AND.  time_stg_call >= dt_stg_call  .AND.                        &
    1245              intermediate_timestep_count == intermediate_timestep_count_max )  THEN
     1247!--       Impose a turbulent inflow using synthetic generated turbulence.
     1248          IF ( use_syn_turb_gen  .AND.                                                             &
     1249               intermediate_timestep_count == intermediate_timestep_count_max )  THEN
    12461250             CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
    12471251             CALL stg_main
Note: See TracChangeset for help on using the changeset viewer.