Ignore:
Timestamp:
Jun 12, 2019 11:52:39 AM (2 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

File:
1 edited

Legend:

Unmodified
Added
Removed
  • 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,                                                           &
Note: See TracChangeset for help on using the changeset viewer.