Ignore:
Timestamp:
Oct 15, 2018 2:21:08 PM (3 years ago)
Author:
suehring
Message:

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

File:
1 edited

Legend:

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

    r3289 r3347  
    2525! -----------------
    2626! $Id$
     27! - Revise structure of init routine
     28! - introduce new parameters to skip STG for some timesteps
     29! - introduce time-dependent parametrization of Reynolds-stress tensor
     30! - Bugfix in allocation of mean_inflow_profiles
     31!
     32! 3341 2018-10-15 10:31:27Z suehring
    2733! Introduce module parameter for number of inflow profiles
    2834!
     
    122128! Authors:
    123129! --------
    124 ! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
     130! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
     131!
    125132!
    126133!
     
    133140!> a time scale for each velocity component. The profiles of length and time
    134141!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
    135 !> tensor are read from file STG_PROFILES.
     142!> tensor can be either read from file STG_PROFILES, or will be parametrized
     143!> within the boundary layer.
    136144!>
    137145!> @todo test restart
     
    149157
    150158    USE arrays_3d,                                                             &
    151         ONLY:  mean_inflow_profiles, u, v, w
     159        ONLY:  mean_inflow_profiles, q, pt, u, v, w, zu, zw
    152160
    153161    USE basic_constants_and_equations_mod,                                     &
    154         ONLY:  pi
     162        ONLY:  g, kappa, pi
    155163
    156164    USE control_parameters,                                                    &
     
    159167
    160168    USE cpulog,                                                                &
    161         ONLY:  cpu_log, log_point, log_point_s
     169        ONLY:  cpu_log, log_point
    162170
    163171    USE indices,                                                               &
     
    170178#endif
    171179
     180    USE nesting_offl_mod,                                                      &
     181        ONLY:  zi_ribulk
     182
    172183    USE pegrid,                                                                &
    173184        ONLY:  comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims
     
    184195
    185196
    186     LOGICAL :: velocity_seed_initialized = .FALSE.  !< true after first call of stg_main
    187     LOGICAL :: use_syn_turb_gen = .FALSE.           !< switch to use synthetic turbulence generator
     197    LOGICAL ::  velocity_seed_initialized = .FALSE.     !< true after first call of stg_main
     198    LOGICAL ::  parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized
     199    LOGICAL ::  use_syn_turb_gen = .FALSE.              !< switch to use synthetic turbulence generator
    188200
    189201    INTEGER(iwp) ::  id_stg_left        !< left lateral boundary core id in case of turbulence generator
     
    202214    INTEGER(iwp) ::  nzt_y_stg          !< upper bound of z coordinate (required for transposing z on PEs along y)
    203215
    204     REAL(wp) :: mc_factor    !< mass flux correction factor
    205 
    206     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz      !< displacement for MPI_GATHERV
    207     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz  !< receive count for MPI_GATHERV
    208     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz      !< displacement for MPI_GATHERV
    209     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz  !< receive count for MPI_GATHERV
    210     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux            !< length scale of u in x direction (in gp)
    211     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy            !< length scale of u in y direction (in gp)
    212     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz            !< length scale of u in z direction (in gp)
    213     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx            !< length scale of v in x direction (in gp)
    214     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy            !< length scale of v in y direction (in gp)
    215     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz            !< length scale of v in z direction (in gp)
    216     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx            !< length scale of w in x direction (in gp)
    217     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy            !< length scale of w in y direction (in gp)
    218     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz            !< length scale of w in z direction (in gp)
    219 
    220     INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed        !< seed of random number for rn-generator
    221 
    222     REAL(wp), DIMENSION(:), ALLOCATABLE :: a11             !< coefficient for Lund rotation
    223     REAL(wp), DIMENSION(:), ALLOCATABLE :: a21             !< coefficient for Lund rotation
    224     REAL(wp), DIMENSION(:), ALLOCATABLE :: a22             !< coefficient for Lund rotation
    225     REAL(wp), DIMENSION(:), ALLOCATABLE :: a31             !< coefficient for Lund rotation
    226     REAL(wp), DIMENSION(:), ALLOCATABLE :: a32             !< coefficient for Lund rotation
    227     REAL(wp), DIMENSION(:), ALLOCATABLE :: a33             !< coefficient for Lund rotation
    228     REAL(wp), DIMENSION(:), ALLOCATABLE :: tu              !< Lagrangian time scale of u
    229     REAL(wp), DIMENSION(:), ALLOCATABLE :: tv              !< Lagrangian time scale of v
    230     REAL(wp), DIMENSION(:), ALLOCATABLE :: tw              !< Lagrangian time scale of w
    231 
    232     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux           !< filter function for u in x direction
    233     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy           !< filter function for u in y direction
    234     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz           !< filter function for u in z direction
    235     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx           !< filter function for v in x direction
    236     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy           !< filter function for v in y direction
    237     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz           !< filter function for v in z direction
    238     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx           !< filter function for w in y direction
    239     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy           !< filter function for w in y direction
    240     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz           !< filter function for w in z direction
    241     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz         !< velocity seed for u at xz plane
    242     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz        !< velocity seed for u at xz plane with new random number
    243     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz         !< velocity seed for u at yz plane
    244     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz        !< velocity seed for u at yz plane with new random number
    245     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz         !< velocity seed for v at xz plane
    246     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz        !< velocity seed for v at xz plane with new random number
    247     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz         !< velocity seed for v at yz plane
    248     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz        !< velocity seed for v at yz plane with new random number
    249     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz         !< velocity seed for w at xz plane
    250     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz        !< velocity seed for w at xz plane with new random number
    251     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz         !< velocity seed for w at yz plane
    252     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz        !< velocity seed for w at yz plane with new random number
    253 
     216    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_xz      !< displacement for MPI_GATHERV
     217    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_xz  !< receive count for MPI_GATHERV
     218    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  displs_yz      !< displacement for MPI_GATHERV
     219    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  recv_count_yz  !< receive count for MPI_GATHERV
     220    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nux            !< length scale of u in x direction (in gp)
     221    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nuy            !< length scale of u in y direction (in gp)
     222    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nuz            !< length scale of u in z direction (in gp)
     223    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvx            !< length scale of v in x direction (in gp)
     224    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvy            !< length scale of v in y direction (in gp)
     225    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nvz            !< length scale of v in z direction (in gp)
     226    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwx            !< length scale of w in x direction (in gp)
     227    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwy            !< length scale of w in y direction (in gp)
     228    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  nwz            !< length scale of w in z direction (in gp)
     229    INTEGER(iwp), DIMENSION(:), ALLOCATABLE ::  seed           !< seed of random number for rn-generator
     230
     231    REAL(wp) ::  cosmo_ref = 10.0_wp      !< height of first vertical grid level in mesoscale model, used for calculation of scaling parameters
     232    REAL(wp) ::  dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics
     233    REAL(wp) ::  dt_stg_call = 5.0_wp     !< time interval for calling synthetic turbulence generator
     234    REAL(wp) ::  mc_factor                !< mass flux correction factor
     235    REAL(wp) ::  scale_l                  !< scaling parameter used for turbulence parametrization - Obukhov length
     236    REAL(wp) ::  scale_us                 !< scaling parameter used for turbulence parametrization - friction velocity
     237    REAL(wp) ::  scale_wm                 !< scaling parameter used for turbulence parametrization - momentum scale 
     238    REAL(wp) ::  time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information   
     239    REAL(wp) ::  time_stg_call = 0.0_wp   !< time counter for calling generator   
     240   
     241   
     242    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r11              !< Reynolds parameter
     243    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r21              !< Reynolds parameter
     244    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r22              !< Reynolds parameter
     245    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r31              !< Reynolds parameter
     246    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r32              !< Reynolds parameter
     247    REAL(wp),DIMENSION(:), ALLOCATABLE ::  r33              !< Reynolds parameter
     248   
     249    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a11             !< coefficient for Lund rotation
     250    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a21             !< coefficient for Lund rotation
     251    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a22             !< coefficient for Lund rotation
     252    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a31             !< coefficient for Lund rotation
     253    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a32             !< coefficient for Lund rotation
     254    REAL(wp), DIMENSION(:), ALLOCATABLE ::  a33             !< coefficient for Lund rotation
     255    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tu              !< Lagrangian time scale of u
     256    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tv              !< Lagrangian time scale of v
     257    REAL(wp), DIMENSION(:), ALLOCATABLE ::  tw              !< Lagrangian time scale of w
     258
     259    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bux           !< filter function for u in x direction
     260    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  buy           !< filter function for u in y direction
     261    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  buz           !< filter function for u in z direction
     262    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvx           !< filter function for v in x direction
     263    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvy           !< filter function for v in y direction
     264    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bvz           !< filter function for v in z direction
     265    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwx           !< filter function for w in y direction
     266    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwy           !< filter function for w in y direction
     267    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  bwz           !< filter function for w in z direction
     268    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fu_xz         !< velocity seed for u at xz plane
     269    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fuo_xz        !< velocity seed for u at xz plane with new random number
     270    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fu_yz         !< velocity seed for u at yz plane
     271    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fuo_yz        !< velocity seed for u at yz plane with new random number
     272    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fv_xz         !< velocity seed for v at xz plane
     273    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fvo_xz        !< velocity seed for v at xz plane with new random number
     274    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fv_yz         !< velocity seed for v at yz plane
     275    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fvo_yz        !< velocity seed for v at yz plane with new random number
     276    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fw_xz         !< velocity seed for w at xz plane
     277    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_xz        !< velocity seed for w at xz plane with new random number
     278    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fw_yz         !< velocity seed for w at yz plane
     279    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_yz        !< velocity seed for w at yz plane with new random number
     280   
     281    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz     !< imposed disturbances at north/south boundary
     282    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz     !< imposed disturbances at north/south boundary
    254283
    255284!
    256285!-- PALM interfaces:
     286!-- Adjust time and lenght scales, Reynolds stress, and filter functions
     287    INTERFACE stg_adjust
     288       MODULE PROCEDURE stg_adjust
     289    END INTERFACE stg_adjust
     290!
    257291!-- Input parameter checks to be done in check_parameters
    258292    INTERFACE stg_check_parameters
     
    319353!
    320354!-- Public interfaces
    321     PUBLIC  stg_check_parameters, stg_header, stg_init, stg_main, stg_parin,   &
    322             stg_wrd_global, stg_rrd_global
     355    PUBLIC  stg_adjust, stg_check_parameters, stg_header, stg_init, stg_main,  &
     356            stg_parin, stg_rrd_global, stg_wrd_global
    323357
    324358!
    325359!-- Public variables
    326     PUBLIC  id_stg_left, id_stg_north, id_stg_right, id_stg_south,             &
    327             use_syn_turb_gen
     360    PUBLIC  dt_stg_call, dt_stg_adjust, id_stg_left, id_stg_north,             &
     361            id_stg_right, id_stg_south, parametrize_inflow_turbulence,         &
     362            time_stg_adjust, time_stg_call, use_syn_turb_gen
    328363
    329364
     
    376411             message_string = 'Using synthetic turbulence generator ' //       &
    377412                              'requires %initializing_actions = '         //   &
    378                               '"set_constant_profiles" or "read_restart_data"'
     413                              '"set_constant_profiles" or "read_restart_data"' //&
     414                              ', if not offline nesting is applied.'
    379415             CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
    380416          ENDIF
     
    382418          IF ( bc_lr /= 'dirichlet/radiation' )  THEN
    383419             message_string = 'Using synthetic turbulence generator ' //       &
    384                               'requires &bc_lr = "dirichlet/radiation"'
     420                              'requires &bc_lr = "dirichlet/radiation", ' //   &
     421                              'if not offline nesting is applied.'
    385422             CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
    386423          ENDIF
    387424          IF ( bc_ns /= 'cyclic' )  THEN
    388425             message_string = 'Using synthetic turbulence generator ' //       &
    389                               'requires &bc_ns = "cyclic"'
     426                              'requires &bc_ns = "cyclic", ' //               &
     427                              'if not offline nesting is applied.'
    390428             CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
    391429          ENDIF
     
    425463       WRITE( io, 3 )
    426464    ENDIF
     465   
     466    IF ( parametrize_inflow_turbulence )  THEN
     467       WRITE( io, 4 ) dt_stg_adjust
     468    ELSE
     469       WRITE( io, 5 )
     470    ENDIF
    427471
    4284721   FORMAT (//' Synthetic turbulence generator information:'/                  &
    429473              ' ------------------------------------------'/)
    430 2   FORMAT ('    synthetic turbulence generator switched on')
    431 3   FORMAT ('    synthetic turbulence generator switched off')
     4742   FORMAT ('    synthetic turbulence generator is switched on')
     4753   FORMAT ('    synthetic turbulence generator is switched off')
     4764   FORMAT ('    imposed turbulence statistics are parametrized and ' \\       &
     477            'ajdusted to boundary-layer development each ', F8.2, ' s' )
     4785   FORMAT ('    imposed turbulence is read from file' )
    432479
    433480 END SUBROUTINE stg_header
     
    443490
    444491    USE arrays_3d,                                                             &
    445         ONLY:  dzw, ddzw, u_init, v_init, zu, zw
     492        ONLY:  dzw, ddzw, u_init, v_init, zu
    446493
    447494    USE control_parameters,                                                    &
     
    460507    IMPLICIT NONE
    461508
    462     LOGICAL ::  file_stg_exist = .FALSE. !< flag indication whether parameter file for Reynolds stress and length scales exist
     509    LOGICAL ::  file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist
    463510
    464511#if defined( __parallel )
     
    484531    REAL(wp) :: d21     !< vertical interpolation for a21
    485532    REAL(wp) :: d22     !< vertical interpolation for a22
    486     REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales
    487533    REAL(wp) :: luy     !< length scale for u in y direction
    488534    REAL(wp) :: luz     !< length scale for u in z direction
     
    494540    REAL(wp) :: zz      !< height
    495541
    496     REAL(wp) :: length_scale_surface !< typical length scale
    497     REAL(wp) :: r_ii_0               !< correction factor
    498     REAL(wp) :: time_scale           !< typical time scale
    499     REAL(wp) :: length_scale_z       !< typical length scale
    500 
    501     REAL(wp),DIMENSION(nzb:nzt+1) :: r11  !< Reynolds parameter
    502     REAL(wp),DIMENSION(nzb:nzt+1) :: r21  !< Reynolds parameter
    503     REAL(wp),DIMENSION(nzb:nzt+1) :: r22  !< Reynolds parameter
    504     REAL(wp),DIMENSION(nzb:nzt+1) :: r31  !< Reynolds parameter
    505     REAL(wp),DIMENSION(nzb:nzt+1) :: r32  !< Reynolds parameter
    506     REAL(wp),DIMENSION(nzb:nzt+1) :: r33  !< Reynolds parameter
    507 
    508542
    509543#if defined( __parallel )
     
    512546
    513547    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
    514 
    515     IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                             &
    516        ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
    517 
    518     ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
    519                a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
    520                nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
    521                nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
    522                nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
    523                tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
    524548
    525549#if defined( __parallel )
     
    615639    ENDDO
    616640    CALL RANDOM_SEED(put=seed)
    617 
     641!
     642!-- Allocate required arrays
     643!-- mean_inflow profiles must not be allocated in offline nesting
     644    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
     645       IF ( .NOT. ALLOCATED( mean_inflow_profiles ) )                          &
     646          ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )
     647    ENDIF
     648
     649    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                 &
     650               a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1),                 &
     651               nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1),                 &
     652               nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1),                 &
     653               nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1),                 &
     654               r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1),                 &
     655               r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1),                 &
     656               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1)   )
     657               
     658    ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) )
     659    ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) )
     660    dist_xz = 0.0_wp
     661    dist_yz = 0.0_wp
     662!
    618663!-- Read inflow profile
    619664!-- Height levels of profiles in input profile is as follows:
     
    669714       
    670715       CLOSE( 90 )
    671 
     716!
     717!--    Calculate coefficient matrix from Reynolds stress tensor 
     718!--    (Lund rotation)
     719       CALL calc_coeff_matrix
     720!
     721!-- No information about turbulence and its length scales are available.
     722!-- Instead, parametrize turbulence which is imposed at the boundaries.
     723!-- Set flag which indicates that turbulence is parametrized, which is done
     724!-- later when energy-balance models are already initialized. This is
     725!-- because the STG needs information about surface properties such as
     726!-- roughness to build 'realistic' turbulence profiles.
    672727    ELSE
    673728!
    674 !--    Set-up default turbulent length scales. From the numerical point of
    675 !--    view the imposed perturbations should not be immediately dissipated
    676 !--    by the numerics. The numerical dissipation, however, acts on scales
    677 !--    up to 8 x the grid spacing. For this reason, set the turbulence
    678 !--    length scale to 8 time the grid spacing. Further, above the boundary
    679 !--    layer height, set turbulence lenght scales to zero (equivalent to not
    680 !--    imposing any perturbations) in order to save computational costs.
    681 !--    Typical length (time) scales of 100 m (s) should be a good compromise
    682 !--    between all stratrifications. Near-surface variances are fixed to
    683 !--    0.1 m2/s2, vertical fluxes are one order of magnitude smaller.
    684 !--    Vertical fluxes
    685        length_scale_surface = 100.0_wp
    686        r_ii_0               = 0.1_wp
    687        time_scale           = 100.0_wp
    688        DO  k = nzb+1, nzt+1
    689           dum_exp        = MERGE( -zu(k) / ( 0.3* zu(nzt) ),                   &
    690                                   0.0_wp,                                      &
    691                                   zu(k) > 0.3 * zu(nzt)                        &
    692                                 )
    693           length_scale_z = 8.0_wp * MIN( dx, dy, dzw(k) )
    694          
    695 !           IF ( zu(k) > zi_richardson ) 
    696 
    697           nux(k) = MAX( INT( length_scale_z * ddx     ), 1 )
    698           nuy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
    699           nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
    700           nvx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
    701           nvy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
    702           nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
    703           nwx(k) = MAX( INT( length_scale_z * ddx     ), 1 )
    704           nwy(k) = MAX( INT( length_scale_z * ddy     ), 1 )
    705           nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 )
    706 
    707           r11(k) = r_ii_0 * EXP( dum_exp )
    708           r22(k) = r_ii_0 * EXP( dum_exp )
    709           r33(k) = r_ii_0 * EXP( dum_exp )
    710 
    711           r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
    712           r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
    713           r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp )
    714 
    715           tu(k)  = time_scale
    716           tv(k)  = time_scale
    717           tw(k)  = time_scale
    718 
    719        ENDDO
    720        
    721        nux(nzb) = nux(nzb+1)
    722        nuy(nzb) = nuy(nzb+1)
    723        nuz(nzb) = nuz(nzb+1)
    724        nvx(nzb) = nvx(nzb+1)
    725        nvy(nzb) = nvy(nzb+1)
    726        nvz(nzb) = nvz(nzb+1)
    727        nwx(nzb) = nwx(nzb+1)
    728        nwy(nzb) = nwy(nzb+1)
    729        nwz(nzb) = nwz(nzb+1)
    730 
    731        r11(nzb) = r11(nzb+1)
    732        r22(nzb) = r22(nzb+1)
    733        r33(nzb) = r33(nzb+1)
    734 
    735        r21(nzb) = r11(nzb+1)
    736        r31(nzb) = r31(nzb+1)
    737        r32(nzb) = r32(nzb+1)
    738 
    739        tu(nzb)  = time_scale
    740        tv(nzb)  = time_scale
    741        tw(nzb)  = time_scale
    742 
     729!--    Set flag indicating that turbulence is parametrized
     730       parametrize_inflow_turbulence = .TRUE.
     731!
     732!--    Determine boundary-layer depth, which is used to initialize lenght
     733!--    scales
     734       CALL calc_scaling_variables
     735!
     736!--    Initialize lenght and time scales, which in turn are used
     737!--    to initialize the filter functions.
     738       CALL calc_length_and_time_scale
     739           
    743740    ENDIF
    744741
    745742!
    746 !-- Assign initial profiles
     743!-- Assign initial profiles. Note, this is only required if turbulent
     744!-- inflow from the left is desired, not in case of any of the
     745!-- nesting (offline or self nesting) approaches.
    747746    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
    748747       u_init = mean_inflow_profiles(:,1)
     
    751750       e_init = MAXVAL( mean_inflow_profiles(:,5) )
    752751    ENDIF
    753 !
    754 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
    755     DO  k = nzb, nzt+1
    756        IF ( r11(k) > 0.0_wp )  THEN
    757           a11(k) = SQRT( r11(k) )
    758           a21(k) = r21(k) / a11(k)
    759        ELSE
    760           a11(k) = 0.0_wp
    761           a21(k) = 0.0_wp
    762        ENDIF
    763 
    764        a22(k) = r22(k) - a21(k)**2
    765        IF ( a22(k) > 0.0_wp )  THEN
    766           a22(k) = SQRT( a22(k) )
    767        ELSE
    768           a22(k) = 0.0_wp
    769        ENDIF
    770 
    771 !
    772 !--    a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
    773 !--    d21, d22) because of different vertical grid
    774        IF ( k .le. nzt )  THEN
    775           d11 = 0.5_wp * ( r11(k) + r11(k+1) )
    776           IF ( d11 > 0.0_wp )  THEN
    777              d11 = SQRT( d11 )
    778              d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
    779              a31(k) = r31(k) / d11
    780           ELSE
    781              d21 = 0.0_wp
    782              a31(k) = 0.0_wp
    783           ENDIF
    784 
    785           d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
    786           IF ( d22 > 0.0_wp )  THEN
    787              a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
    788           ELSE
    789              a32(k) = 0.0_wp
    790           ENDIF
    791 
    792           a33(k) = r33(k) - a31(k) ** 2 - a32(k) ** 2
    793           IF ( a33(k) > 0.0_wp )  THEN
    794              a33(k) = SQRT( a33(k) )
    795           ELSE
    796              a33(k) = 0.0_wp
    797           ENDIF
    798        ELSE
    799           a31(k) = a31(k-1)
    800           a32(k) = a32(k-1)
    801           a33(k) = a33(k-1)
    802        ENDIF
    803 
    804     ENDDO
     752   
    805753!
    806754!-- Define the size of the filter functions and allocate them.
     
    961909    ENDDO
    962910
    963 END SUBROUTINE stg_filter_func
     911 END SUBROUTINE stg_filter_func
    964912
    965913
     
    977925
    978926
    979     NAMELIST /stg_par/  use_syn_turb_gen
     927    NAMELIST /stg_par/  dt_stg_adjust, dt_stg_call, use_syn_turb_gen
    980928
    981929    line = ' '
     
    11081056    REAL(wp) :: volume_flow_l   !< local mass flux through lateral boundary
    11091057
    1110     REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary
    1111     REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary
    1112 
    1113 
    11141058    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
    11151059
    11161060!
    11171061!-- Calculate time step which is needed for filter functions
    1118     dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
    1119 
     1062    dt_stg = MAX( dt_3d, dt_stg_call ) !* weight_substep(intermediate_timestep_count)
    11201063!
    11211064!-- Initial value of fu, fv, fw
     
    11451088       velocity_seed_initialized = .TRUE.
    11461089    ENDIF
     1090   
    11471091!
    11481092!-- New set of fu, fv, fw
     
    11501094    CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left )
    11511095    CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left )
    1152 
     1096   
    11531097    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent         &
    1154                              .AND.  .NOT.  rans_mode ) )  THEN
    1155 !
    1156 !--       Generate turbulence at right boundary
    1157           CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
    1158           CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
    1159           CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
    1160 !
    1161 !--       Generate turbulence at north boundary
    1162           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
    1163           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
    1164           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
    1165 !
    1166 !--       Generate turbulence at south boundary
    1167           CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
    1168           CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
    1169           CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
     1098                          .AND.  .NOT.  rans_mode ) )  THEN
     1099!
     1100!--    Generate turbulence at right boundary
     1101       CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right )
     1102       CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right )
     1103       CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right )
     1104!
     1105!--    Generate turbulence at north boundary
     1106       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north )
     1107       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north )
     1108       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north )
     1109!
     1110!--    Generate turbulence at south boundary
     1111       CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south )
     1112       CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south )
     1113       CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south )
    11701114    ENDIF
     1115   
    11711116!
    11721117!-- Turbulence generation at left and or right boundary
    11731118    IF ( myidx == id_stg_left  .OR.  myidx == id_stg_right )  THEN
    1174 
     1119   
    11751120       DO  j = nysg, nyng
    11761121          DO  k = nzb, nzt + 1
     
    11801125                fu_yz(k,j) = fuo_yz(k,j)
    11811126             ELSE
    1182                 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) +     &
     1127                fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
    11831128                         fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
    11841129             ENDIF
     
    11871132                fv_yz(k,j) = fvo_yz(k,j)
    11881133             ELSE
    1189                 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
     1134                fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
    11901135                         fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    11911136             ENDIF
     
    11941139                fw_yz(k,j) = fwo_yz(k,j)
    11951140             ELSE
    1196                 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
     1141                fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
    11971142                         fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    11981143             ENDIF
     
    12041149                dist_yz(k,j,2) = 0.0_wp
    12051150                dist_yz(k,j,3) = 0.0_wp
    1206 !                 dist_yz(k,j,4) = 0.0_wp
    1207 !                 dist_yz(k,j,5) = 0.0_wp
    12081151             ELSE
    1209                 dist_yz(k,j,1) = a11(k) * fu_yz(k,j)
     1152                dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp )
    12101153                !experimental test of 1.2
    1211                 dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
    1212                                          * 1.2_wp )                            &
    1213                                        * (   a21(k) * fu_yz(k,j)               &
    1214                                            + a22(k) * fv_yz(k,j) )
    1215                 dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
    1216                                          * 1.3_wp )                            &
    1217                                        * (   a31(k) * fu_yz(k,j)               &
    1218                                            + a32(k) * fv_yz(k,j)               &
    1219                                            + a33(k) * fw_yz(k,j) )
    1220                 ! Calculation for pt and e not yet implemented
    1221 !                 dist_yz(k,j,4) = 0.0_wp
    1222 !                 dist_yz(k,j,5) = 0.0_wp
     1154                dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
     1155                                      * 1.2_wp )                               &
     1156                                     * (   a21(k) * fu_yz(k,j)                 &
     1157                                         + a22(k) * fv_yz(k,j) ), 3.0_wp )
     1158                dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
     1159                                      * 1.3_wp )                               &
     1160                                    * (   a31(k) * fu_yz(k,j)                  &
     1161                                        + a32(k) * fv_yz(k,j)                  &
     1162                                        + a33(k) * fw_yz(k,j) ), 3.0_wp )
    12231163             ENDIF
    12241164
     
    12411181
    12421182#if defined( __parallel )
    1243           CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
    1244                               1, MPI_REAL, MPI_SUM, comm1dy, ierr )
     1183          CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM,    &
     1184                              comm1dy, ierr )
    12451185#else
    12461186          mc_factor = mc_factor_l
     
    12981238!--       Add disturbances
    12991239          IF ( myidx == id_stg_left  )  THEN
    1300 
    13011240             DO  j = nys, nyn
    1302                 DO  k = nzb, nzt
     1241                DO  k = nzb+1, nzt
    13031242                   u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) )                    &
    13041243                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
    13051244                                            BTEST( wall_flags_0(k,j,0), 1 ) )
     1245                   u(k,j,-1) = u(k,j,0)
    13061246                   v(k,j,-1)  = ( v(k,j,-1)  + dist_yz(k,j,2)  )               &
    13071247                       * mc_factor * MERGE( 1.0_wp, 0.0_wp,                    &
     
    13141254          ENDIF
    13151255          IF ( myidx == id_stg_right  )  THEN
    1316 
    13171256             DO  j = nys, nyn
    1318                 DO  k = nzb, nzt
     1257                DO  k = nzb+1, nzt
    13191258                   u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) )            &
    13201259                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
     
    13351274!-- Turbulence generation at north and south boundary
    13361275    IF ( myidy == id_stg_north  .OR.  myidy == id_stg_south )  THEN
    1337 
     1276   
    13381277       DO  i = nxlg, nxrg
    13391278          DO  k = nzb, nzt + 1
     
    13511290             ELSE
    13521291                fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) +     &
    1353                          fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
     1292                      fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
    13541293             ENDIF
    1355 
     1294           
    13561295             IF ( tw(k) == 0.0_wp )  THEN
    13571296                fw_xz(k,i) = fwo_xz(k,i)
    13581297             ELSE
    13591298                fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) +     &
    1360                          fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
     1299                      fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
    13611300             ENDIF
    13621301!
     
    13691308
    13701309             ELSE
    1371                 dist_xz(k,i,1) = a11(k) * fu_xz(k,i)
     1310                dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp )
    13721311                !experimental test of 1.2
    1373                 dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )                &
    1374                                          * 1.2_wp )                            &
    1375                                        * (   a21(k) * fu_xz(k,i)               &
    1376                                            + a22(k) * fv_xz(k,i) )
    1377                 dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) )                 &
    1378                                          * 1.3_wp )                            &
    1379                                        * (   a31(k) * fu_xz(k,i)               &
    1380                                            + a32(k) * fv_xz(k,i)               &
    1381                                            + a33(k) * fw_xz(k,i) )
     1312                dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) )           &
     1313                                      * 1.2_wp )                               &
     1314                                     * (   a21(k) * fu_xz(k,i)                 &
     1315                                         + a22(k) * fv_xz(k,i) ), 3.0_wp )
     1316                dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) )            &
     1317                                      * 1.3_wp )                               &
     1318                                    * (   a31(k) * fu_xz(k,i)                  &
     1319                                        + a32(k) * fv_xz(k,i)                  &
     1320                                        + a33(k) * fw_xz(k,i) ), 3.0_wp )
    13821321             ENDIF
    13831322
    13841323          ENDDO
    13851324       ENDDO
     1325
    13861326!
    13871327!--    Mass flux correction following Kim et al. (2013)
     
    14251365
    14261366          DO  i = nxl, nxr
    1427              DO  k = nzb, nzt
     1367             DO  k = nzb+1, nzt
    14281368                u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) )                     &
    14291369                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
     
    14321372                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
    14331373                                           BTEST( wall_flags_0(k,0,i), 2 ) )
     1374                v(k,-1,i) = v(k,0,i)
    14341375                w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3)  )                    &
    14351376                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                     &
     
    14411382
    14421383          DO  i = nxl, nxr
    1443              DO  k = nzb, nzt
    1444                 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )                &
    1445                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
     1384             DO  k = nzb+1, nzt
     1385                u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) )               &
     1386                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    14461387                                          BTEST( wall_flags_0(k,nyn+1,i), 1 ) )
    1447                 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,1) )                &
    1448                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
     1388                v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) )               &
     1389                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    14491390                                          BTEST( wall_flags_0(k,nyn+1,i), 2 ) )
    1450                 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,1) )                &
    1451                      * mc_factor * MERGE( 1.0_wp, 0.0_wp,                       &
     1391                w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) )               &
     1392                     * mc_factor * MERGE( 1.0_wp, 0.0_wp,                      &
    14521393                                          BTEST( wall_flags_0(k,nyn+1,i), 3 ) )
    14531394             ENDDO
     
    14551396       ENDIF
    14561397    ENDIF
     1398!
     1399!-- Finally, set time counter for calling STG to zero
     1400    time_stg_call = 0.0_wp
    14571401
    14581402    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
     
    14931437    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
    14941438
    1495     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter func in y-dir
    1496     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
     1439    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_y     !< filter function in y-direction
     1440    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
    14971441    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l   !<  local velocity seed
    14981442    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng)     :: f_n     !<  velocity seed
     
    16131557    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
    16141558    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
    1615     INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_y*2
     1559    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2   !< n_x*2
    16161560    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2   !< n_z*2
    16171561
     
    16201564    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
    16211565
    1622     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter func in y-dir
    1623     REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter func in z-dir
     1566    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_x     !< filter function in x-direction
     1567    REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1)    :: b_z     !< filter function in z-direction
    16241568    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l   !<  local velocity seed
    16251569    REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg)     :: f_n     !<  velocity seed
     
    17121656 END SUBROUTINE stg_generate_seed_xz
    17131657
     1658!------------------------------------------------------------------------------!
     1659! Description:
     1660! ------------
     1661!> Parametrization of the Reynolds stress tensor, following the parametrization
     1662!> described in Rotach et al. (1996), which is applied in state-of-the-art
     1663!> dispserion modelling. Please note, the parametrization does not distinguish
     1664!> between along-wind and cross-wind turbulence.
     1665!------------------------------------------------------------------------------!
     1666 SUBROUTINE parametrize_reynolds_stress
     1667
     1668    USE arrays_3d,                                                             &
     1669        ONLY:  zu
     1670
     1671    IMPLICIT NONE
     1672
     1673    INTEGER(iwp) :: k            !< loop index in z-direction
     1674   
     1675    REAL(wp)     ::  zzi         !< ratio of z/zi
     1676   
     1677!
     1678!--
     1679    DO  k = nzb+1, nzt+1
     1680         
     1681       IF ( zu(k) <= zi_ribulk )  THEN
     1682!
     1683!--       Determine normalized height coordinate
     1684          zzi = zu(k) / zi_ribulk
     1685!
     1686!--       u'u' and v'v'. Assume isotropy. Note, add a small negative number
     1687!--       to the denominator, else the merge-function can crash if scale_l is
     1688!--       zero.
     1689          r11(k) = scale_us**2 * (                                             &
     1690                      MERGE( 0.35_wp * (                                       &
     1691                           - zi_ribulk / ( kappa * scale_l - 10E-4_wp )        &
     1692                                       )**( 2.0_wp / 3.0_wp ),                 &
     1693                             0.0_wp,                                           &
     1694                             scale_l < 0.0_wp )                                &
     1695                    + 5.0_wp - 4.0_wp * zzi                                    &
     1696                                 )
     1697                                 
     1698          r22(k) = r11(k)
     1699!
     1700!--       w'w'
     1701          r33(k) = scale_wm**2 * (                                             &
     1702                      1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) &
     1703                    + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2            &                     
     1704                                 )
     1705!
     1706!--       u'w' and v'w'. Assume isotropy.
     1707          r31(k) = - scale_us**2 * (                                           &
     1708                         1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) )             &
     1709                                   )
     1710                             
     1711          r32(k) = r31(k)
     1712!
     1713!--       For u'v' no parametrization exist so far - ?. For simplicity assume
     1714!--       a similar profile as for u'w'.
     1715          r21(k) = r31(k)
     1716!
     1717!--    Above the boundary layer, assmume laminar flow conditions.
     1718       ELSE
     1719          r11(k) = 10E-8_wp
     1720          r22(k) = 10E-8_wp
     1721          r33(k) = 10E-8_wp
     1722          r21(k) = 10E-8_wp
     1723          r31(k) = 10E-8_wp
     1724          r32(k) = 10E-8_wp
     1725       ENDIF
     1726!        write(9,*) zu(k), r11(k), r33(k), r31(k), zi_ribulk, scale_us, scale_wm, scale_l
     1727    ENDDO
     1728
     1729!
     1730!-- Set bottom boundary condition   
     1731    r11(nzb) = r11(nzb+1)
     1732    r22(nzb) = r22(nzb+1)
     1733    r33(nzb) = r33(nzb+1)
     1734
     1735    r21(nzb) = r11(nzb+1)
     1736    r31(nzb) = r31(nzb+1)
     1737    r32(nzb) = r32(nzb+1)
     1738   
     1739
     1740 END SUBROUTINE parametrize_reynolds_stress
     1741 
     1742!------------------------------------------------------------------------------!
     1743! Description:
     1744! ------------
     1745!> Calculate the coefficient matrix from the Lund rotation.
     1746!------------------------------------------------------------------------------!
     1747 SUBROUTINE calc_coeff_matrix
     1748
     1749    IMPLICIT NONE
     1750
     1751    INTEGER(iwp) :: k   !< loop index in z-direction
     1752   
     1753!
     1754!-- Calculate coefficient matrix. Split loops to allow for loop vectorization.
     1755    DO  k = nzb+1, nzt+1
     1756       IF ( r11(k) > 0.0_wp )  THEN
     1757          a11(k) = SQRT( r11(k) )
     1758          a21(k) = r21(k) / a11(k)
     1759          a31(k) = r31(k) / a11(k)
     1760       ELSE
     1761          a11(k) = 10E-8_wp
     1762          a21(k) = 10E-8_wp
     1763          a31(k) = 10E-8_wp
     1764       ENDIF
     1765    ENDDO
     1766    DO  k = nzb+1, nzt+1
     1767       a22(k) = r22(k) - a21(k)**2
     1768       IF ( a22(k) > 0.0_wp )  THEN
     1769          a22(k) = SQRT( a22(k) )
     1770          a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
     1771       ELSE
     1772          a22(k) = 10E-8_wp
     1773          a32(k) = 10E-8_wp
     1774       ENDIF
     1775    ENDDO
     1776    DO  k = nzb+1, nzt+1
     1777       a33(k) = r33(k) - a31(k)**2 - a32(k)**2
     1778       IF ( a33(k) > 0.0_wp )  THEN
     1779          a33(k) =  SQRT( a33(k) )
     1780       ELSE
     1781          a33(k) = 10E-8_wp
     1782       ENDIF
     1783    ENDDO   
     1784!
     1785!-- Set bottom boundary condition
     1786    a11(nzb) = a11(nzb+1)
     1787    a22(nzb) = a22(nzb+1)
     1788    a21(nzb) = a21(nzb+1)
     1789    a33(nzb) = a33(nzb+1)   
     1790    a31(nzb) = a31(nzb+1)
     1791    a32(nzb) = a32(nzb+1)   
     1792
     1793 END SUBROUTINE calc_coeff_matrix
     1794 
     1795!------------------------------------------------------------------------------!
     1796! Description:
     1797! ------------
     1798!> This routine controls the re-adjustment of the turbulence statistics used
     1799!> for generating turbulence at the lateral boundaries. 
     1800!------------------------------------------------------------------------------!
     1801 SUBROUTINE stg_adjust
     1802
     1803    IMPLICIT NONE
     1804   
     1805    INTEGER(iwp) ::  k
     1806
     1807    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
     1808!
     1809!-- Compute mean boundary layer height according to Richardson-Bulk
     1810!-- criterion using the inflow profiles. Further velocity scale as well as
     1811!-- mean friction velocity are calculated.
     1812    CALL calc_scaling_variables
     1813!
     1814!-- Set length and time scales depending on boundary-layer height
     1815    CALL calc_length_and_time_scale
     1816!
     1817!-- Parametrize Reynolds-stress tensor, diagonal elements as well
     1818!-- as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows
     1819!-- Rotach et al. (1996) and is based on boundary-layer depth,
     1820!-- friction velocity and velocity scale.
     1821    CALL parametrize_reynolds_stress
     1822!
     1823!-- Calculate coefficient matrix from Reynolds stress tensor 
     1824!-- (Lund rotation)
     1825    CALL calc_coeff_matrix
     1826!
     1827!-- Determine filter functions on basis of updated length scales
     1828    CALL stg_filter_func( nux, bux ) !filter ux
     1829    CALL stg_filter_func( nuy, buy ) !filter uy
     1830    CALL stg_filter_func( nuz, buz ) !filter uz
     1831    CALL stg_filter_func( nvx, bvx ) !filter vx
     1832    CALL stg_filter_func( nvy, bvy ) !filter vy
     1833    CALL stg_filter_func( nvz, bvz ) !filter vz
     1834    CALL stg_filter_func( nwx, bwx ) !filter wx
     1835    CALL stg_filter_func( nwy, bwy ) !filter wy
     1836    CALL stg_filter_func( nwz, bwz ) !filter wz
     1837!
     1838!-- Reset time counter for controlling next adjustment to zero
     1839    time_stg_adjust = 0.0_wp
     1840   
     1841    CALL  cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
     1842   
     1843 END SUBROUTINE stg_adjust
     1844 
     1845 
     1846!------------------------------------------------------------------------------!
     1847! Description:
     1848! ------------
     1849!> Calculates turbuluent length and time scales if these are not available
     1850!> from measurements.
     1851!------------------------------------------------------------------------------!
     1852 SUBROUTINE calc_length_and_time_scale
     1853
     1854    USE arrays_3d,                                                             &
     1855        ONLY:  dzw, ddzw, u_init, v_init, zu, zw
     1856
     1857    USE grid_variables,                                                        &
     1858        ONLY:  ddx, ddy, dx, dy
     1859       
     1860    USE statistics,                                                            &
     1861        ONLY:  hom
     1862
     1863    IMPLICIT NONE
     1864
     1865
     1866    INTEGER(iwp) :: k            !< loop index in z-direction
     1867
     1868    REAL(wp) ::  length_scale    !< typical length scale
     1869   
     1870!
     1871!-- In initial call the boundary-layer depth can be zero. This case, set
     1872!-- minimum value for boundary-layer depth, to setup length scales correctly.
     1873    zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
     1874!
     1875!-- Set-up default turbulent length scales. From the numerical point of
     1876!-- view the imposed perturbations should not be immediately dissipated
     1877!-- by the numerics. The numerical dissipation, however, acts on scales
     1878!-- up to 8 x the grid spacing. For this reason, set the turbulence
     1879!-- length scale to 8 time the grid spacing. Further, above the boundary
     1880!-- layer height, set turbulence lenght scales to zero (equivalent to not
     1881!-- imposing any perturbations) in order to save computational costs.
     1882!-- Typical time scales are derived by assuming Taylors's hypothesis,
     1883!-- using the length scales and the mean profiles of u- and v-component.
     1884    DO  k = nzb+1, nzt+1
     1885   
     1886       length_scale = 8.0_wp * MIN( dx, dy, dzw(k) )
     1887         
     1888       IF ( zu(k) <= zi_ribulk )  THEN 
     1889!
     1890!--       Assume isotropic turbulence length scales
     1891          nux(k) = MAX( INT( length_scale * ddx     ), 1 )
     1892          nuy(k) = MAX( INT( length_scale * ddy     ), 1 )
     1893          nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
     1894          nvx(k) = MAX( INT( length_scale * ddx     ), 1 )
     1895          nvy(k) = MAX( INT( length_scale * ddy     ), 1 )
     1896          nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
     1897          nwx(k) = MAX( INT( length_scale * ddx     ), 1 )
     1898          nwy(k) = MAX( INT( length_scale * ddy     ), 1 )
     1899          nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 )
     1900!
     1901!--       Limit time scales, else they become very larger for low wind speed,
     1902!--       imposing long-living inflow perturbations which in turn propagate
     1903!--       further into the model domain. Use u_init and v_init to calculate
     1904!--       the time scales, which will be equal to the inflow profiles, both,
     1905!--       in offline nesting mode or in dirichlet/radiation mode.
     1906          tu(k)  = MIN( dt_stg_adjust, length_scale /                          &
     1907                                  ( ABS( u_init(k) ) + 0.1_wp ) )
     1908          tv(k)  = MIN( dt_stg_adjust, length_scale /                          &
     1909                                  ( ABS( v_init(k) ) + 0.1_wp ) )
     1910!
     1911!--       Time scale of w-component is a mixture from u- and v-component.
     1912          tw(k)  = SQRT( tu(k)**2 + tv(k)**2 )
     1913!
     1914!--    Above the boundary layer length scales are zero, i.e. imposed turbulence
     1915!--    is not correlated in space and time, just white noise. This saves
     1916!--    computations power.
     1917       ELSE
     1918          nux(k) = 0.0_wp
     1919          nuy(k) = 0.0_wp
     1920          nuz(k) = 0.0_wp
     1921          nvx(k) = 0.0_wp
     1922          nvy(k) = 0.0_wp
     1923          nvz(k) = 0.0_wp
     1924          nwx(k) = 0.0_wp
     1925          nwy(k) = 0.0_wp
     1926          nwz(k) = 0.0_wp
     1927         
     1928          tu(k)  = 0.0_wp
     1929          tv(k)  = 0.0_wp
     1930          tw(k)  = 0.0_wp
     1931       ENDIF
     1932    ENDDO
     1933!
     1934!-- Set bottom boundary condition for the length and time scales
     1935    nux(nzb) = nux(nzb+1)
     1936    nuy(nzb) = nuy(nzb+1)
     1937    nuz(nzb) = nuz(nzb+1)
     1938    nvx(nzb) = nvx(nzb+1)
     1939    nvy(nzb) = nvy(nzb+1)
     1940    nvz(nzb) = nvz(nzb+1)
     1941    nwx(nzb) = nwx(nzb+1)
     1942    nwy(nzb) = nwy(nzb+1)
     1943    nwz(nzb) = nwz(nzb+1)
     1944
     1945    tu(nzb)  = tu(nzb+1)
     1946    tv(nzb)  = tv(nzb+1)
     1947    tw(nzb)  = tw(nzb+1)
     1948
     1949
     1950 END SUBROUTINE calc_length_and_time_scale
     1951
     1952!------------------------------------------------------------------------------!
     1953! Description:
     1954! ------------
     1955!> Calculate scaling variables which are used for turbulence parametrization
     1956!> according to Rotach et al. (1996). Scaling variables are: friction velocity,
     1957!> boundary-layer depth, momentum velocity scale, and Obukhov length.
     1958!------------------------------------------------------------------------------!
     1959 SUBROUTINE calc_scaling_variables
     1960
     1961
     1962    USE control_parameters,                                                    &
     1963        ONLY:  bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, &
     1964               humidity, pt_surface
     1965             
     1966    USE indices,                                                               &
     1967        ONLY:  nx, ny 
     1968       
     1969    USE surface_mod,                                                           &
     1970        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
     1971
     1972    IMPLICIT NONE
     1973
     1974    INTEGER(iwp) :: i            !< loop index in x-direction
     1975    INTEGER(iwp) :: j            !< loop index in y-direction
     1976    INTEGER(iwp) :: k            !< loop index in z-direction
     1977    INTEGER(iwp) :: k_ref        !< index in z-direction for reference height
     1978    INTEGER(iwp) :: m            !< surface element index
     1979
     1980    REAL(wp) ::  friction_vel_l         !< mean friction veloctiy on subdomain
     1981    REAL(wp) ::  shf_mean               !< mean surface sensible heat flux
     1982    REAL(wp) ::  shf_mean_l             !< mean surface sensible heat flux on subdomain
     1983    REAL(wp) ::  u_int                  !< u-component
     1984    REAL(wp) ::  v_int                  !< v-component
     1985    REAL(wp) ::  vpt_surface            !< near-surface virtual potential temperature
     1986    REAL(wp) ::  w_convective           !< convective velocity scale
     1987    REAL(wp) ::  z0_mean                !< mean roughness length
     1988    REAL(wp) ::  z0_mean_l              !< mean roughness length on subdomain
     1989   
     1990!
     1991!-- Mean friction velocity and velocity scale. Therefore,
     1992!-- pre-calculate mean roughness length and surface sensible heat flux
     1993!-- in the model domain, which are further used to estimate friction
     1994!-- velocity and velocity scale. Note, for z0 linear averaging is applied,
     1995!-- even though this is known to unestimate the effective roughness.
     1996!-- This need to be revised later.
     1997    z0_mean_l  = 0.0_wp
     1998    shf_mean_l = 0.0_wp
     1999    DO  m = 1, surf_def_h(0)%ns
     2000       z0_mean_l  = z0_mean_l  + surf_def_h(0)%z0(m)
     2001       shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m)
     2002    ENDDO   
     2003    DO  m = 1, surf_lsm_h%ns
     2004       z0_mean_l  = z0_mean_l  + surf_lsm_h%z0(m)
     2005       shf_mean_l = shf_mean_l + surf_lsm_h%shf(m)
     2006    ENDDO
     2007    DO  m = 1, surf_usm_h%ns
     2008       z0_mean_l  = z0_mean_l  + surf_usm_h%z0(m)
     2009       shf_mean_l = shf_mean_l + surf_usm_h%shf(m)
     2010    ENDDO
     2011   
     2012#if defined( __parallel )
     2013    CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM,              &
     2014                        comm2d, ierr )
     2015    CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM,            &
     2016                        comm2d, ierr )
     2017#else
     2018    z0_mean  = z0_mean_l
     2019    shf_mean = shf_mean_l
     2020#endif
     2021    z0_mean  = z0_mean  / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
     2022    shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
     2023       
     2024!
     2025!-- Note, Inifor does not use logarithmic interpolation of the
     2026!-- velocity components near the ground, so that near-surface
     2027!-- wind speed and thus the friction velocity is overestimated.
     2028!-- However, friction velocity is used for turbulence
     2029!-- parametrization, so that more physically meaningful values are important.
     2030!-- Hence, derive friction velocity from wind speed at a reference height.
     2031!-- For a first guess use 20 m, which is in the range of the first
     2032!-- COSMO vertical level.
     2033    k_ref = MINLOC( ABS( zu - cosmo_ref ), DIM = 1 ) - 1
     2034   
     2035    friction_vel_l = 0.0_wp
     2036    IF ( bc_dirichlet_l  .OR.  bc_dirichlet_r )  THEN
     2037       
     2038       i = MERGE( -1, nxr + 1, bc_dirichlet_l )
     2039
     2040       DO  j = nys, nyn
     2041!
     2042!--       Note, in u- and v- component the imposed perturbations
     2043!--       from the STG are already included. Check whether this
     2044!--       makes any difference compared to using the pure-mean
     2045!--       inflow profiles.
     2046          u_int = MERGE( u(k_ref,j,i+1), u(k_ref,j,i), bc_dirichlet_l )
     2047          v_int = v(k_ref,j,i)
     2048!
     2049!--       Calculate friction velocity and sum-up. Therefore, assume
     2050!--       neutral condtions.
     2051          friction_vel_l = friction_vel_l + kappa *                            &
     2052                            SQRT( u_int * u_int + v_int * v_int )  /           &
     2053                            LOG( zu(k_ref) / z0_mean )
     2054         
     2055       ENDDO
     2056       
     2057    ENDIF
     2058   
     2059    IF ( bc_dirichlet_s  .OR.  bc_dirichlet_n )  THEN
     2060       
     2061       j = MERGE( -1, nyn + 1, bc_dirichlet_s )
     2062   
     2063       DO  i = nxl, nxr
     2064         
     2065          u_int = u(k_ref,j,i)
     2066          v_int = MERGE( v(k_ref,j+1,i), v(k_ref,j,i), bc_dirichlet_s )
     2067 
     2068          friction_vel_l = friction_vel_l + kappa *                            &
     2069                            SQRT( u_int * u_int + v_int * v_int )  /           &
     2070                            LOG( zu(k_ref) / z0_mean )                       
     2071         
     2072       ENDDO
     2073       
     2074    ENDIF
     2075   
     2076#if defined( __parallel )
     2077    CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM,        &
     2078                        comm2d, ierr )
     2079#else
     2080    scale_us = friction_vel_l
     2081#endif
     2082    scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp )
     2083   
     2084!
     2085!-- Compute mean Obukhov length
     2086    IF ( shf_mean > 0.0_wp )  THEN
     2087       scale_l = - pt_surface / ( kappa * g ) * scale_us**3 / shf_mean
     2088    ELSE
     2089       scale_l = 0.0_wp
     2090    ENDIF
     2091   
     2092!
     2093!-- Compute mean convective velocity scale. Note, in case the mean heat flux
     2094!-- is negative, set convective velocity scale to zero.
     2095    IF ( shf_mean > 0.0_wp )  THEN
     2096       w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp )
     2097    ELSE
     2098       w_convective = 0.0_wp
     2099    ENDIF
     2100!
     2101!-- Finally, in order to consider also neutral or stable stratification,
     2102!-- compute momentum velocity scale from u* and convective velocity scale,
     2103!-- according to Rotach et al. (1996).
     2104    scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp )
     2105   
     2106 END SUBROUTINE calc_scaling_variables
     2107
    17142108 END MODULE synthetic_turbulence_generator_mod
Note: See TracChangeset for help on using the changeset viewer.