Changeset 4566 for palm


Ignore:
Timestamp:
Jun 16, 2020 10:11:51 AM (6 months ago)
Author:
suehring
Message:

synthetic turbulence generator: revise parametrization for reynolds-stress components, turbulent length- and time scales; revise computation of velocity disturbances to be consistent to Xie and Castro (2008); change default value of time interval to adjust turbulence parametrization; bugfix in computation of amplitude-tensor (vertical flux of horizontal momentum)

File:
1 edited

Legend:

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

    r4562 r4566  
    2525! -----------------
    2626! $Id$
     27! - revise parametrization for reynolds-stress components, turbulent length- and time scales
     28! - revise computation of velocity disturbances to be consistent to Xie and Castro (2008)
     29! - change default value of time interval to adjust turbulence parametrization
     30! - bugfix in computation of amplitude-tensor (vertical flux of horizontal momentum)
     31!
     32! 4562 2020-06-12 08:38:47Z raasch
    2733! Parts of r4559 re-formatted
    2834!
     
    151157!> parametrized within the boundary layer.
    152158!>
    153 !> @todo test restart
    154 !>       Enable cyclic_fill
     159!> @todo Enable cyclic_fill
    155160!>       Implement turbulence generation for e and pt
    156 !> @todo Input of height-constant length scales via namelist
    157161!> @note <Enter notes on the module>
    158162!> @bug  Height information from input file is not used. Profiles from input must match with current
     
    287291    INTEGER(iwp) ::  id_stg_right       !< right lateral boundary core id in case of turbulence generator
    288292    INTEGER(iwp) ::  id_stg_south       !< south lateral boundary core id in case of turbulence generator
    289     INTEGER(iwp) ::  mergp              !< maximum length scale (in gp)
     293    INTEGER(iwp) ::  k_zi               !< vertical index of boundary-layer top
     294    INTEGER(iwp) ::  mergp_limit = 1000 !< maximum length scale (in gp)
     295    INTEGER(iwp) ::  mergp_x            !< maximum length scale in x (in gp)
     296    INTEGER(iwp) ::  mergp_xy           !< maximum horizontal length scale (in gp)
     297    INTEGER(iwp) ::  mergp_y            !< maximum length scale in y (in gp)
     298    INTEGER(iwp) ::  mergp_z            !< maximum length scale in z (in gp)
    290299    INTEGER(iwp) ::  nzb_x_stg          !< lower bound of z coordinate (required for transposing z on PEs along x)
    291300    INTEGER(iwp) ::  nzt_x_stg          !< upper bound of z coordinate (required for transposing z on PEs along x)
     
    326335
    327336
     337    LOGICAL ::  adjustment_step               = .FALSE.  !< control flag indicating that time and lenght scales have been updated and
     338                                                         !< no time correlation to the timestep before should be considered
    328339    LOGICAL ::  compute_velocity_seeds_local  = .TRUE.   !< switch to decide whether velocity seeds are computed locally
    329340                                                         !< or if computation is distributed over several processes
     
    336347    REAL(wp) ::  blend                     !< value to create gradually and smooth blending of Reynolds stress and length
    337348                                           !< scales above the boundary layer
    338     REAL(wp) ::  blend_coeff = -2.3_wp     !< coefficient used to ensure that blending functions decreases to 1/10 after
     349    REAL(wp) ::  blend_coeff = -9.3_wp     !< coefficient used to ensure that blending functions decreases to 1/10 after
    339350                                           !< one length scale above ABL top
    340351    REAL(wp) ::  d_l                       !< blend_coeff/length_scale
    341     REAL(wp) ::  dt_stg_adjust = 300.0_wp  !< time interval for adjusting turbulence statistics
     352    REAL(wp) ::  d_nxy                     !< inverse of the total number of xy-grid points
     353    REAL(wp) ::  dt_stg_adjust = 1800.0_wp !< time interval for adjusting turbulence statistics
    342354    REAL(wp) ::  dt_stg_call = 0.0_wp      !< time interval for calling synthetic turbulence generator
    343     REAL(wp) ::  length_scale              !< length scale, default is 8 x minimum grid spacing
    344355    REAL(wp) ::  scale_l                   !< scaling parameter used for turbulence parametrization - Obukhov length
    345356    REAL(wp) ::  scale_us                  !< scaling parameter used for turbulence parametrization - friction velocity
     
    390401    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  fwo_yz  !< velocity seed for w at yz plane with new random number
    391402
    392     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz  !< imposed disturbances at north/south boundary
    393     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz  !< imposed disturbances at north/south boundary
     403    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_xz !< disturbances for parallel/crosswind/vertical component at north/south boundary
     404    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  dist_yz !< disturbances for parallel/crosswind/vertical component at north/south boundary
    394405
    395406!
     
    756767       ENDIF
    757768    ENDIF
     769!
     770!-- Assign initial profiles. Note, this is only required if turbulent inflow from the left is
     771!-- desired, not in case of any of the nesting (offline or self nesting) approaches.
     772    IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
     773       u_init = mean_inflow_profiles(:,1)
     774       v_init = mean_inflow_profiles(:,2)
     775       e_init = MAXVAL( mean_inflow_profiles(:,5) )
     776    ENDIF
    758777
    759778    ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1),                                     &
     
    766785               tu(nzb:nzt+1),  tv(nzb:nzt+1),  tw(nzb:nzt+1) )
    767786
     787    r11 = 0.0_wp
     788    r21 = 0.0_wp
     789    r22 = 0.0_wp
     790    r31 = 0.0_wp
     791    r32 = 0.0_wp
     792    r33 = 0.0_wp
     793    tu  = 0.0_wp
     794    tv  = 0.0_wp
     795    tw  = 0.0_wp
     796
    768797    ALLOCATE ( dist_xz(nzb:nzt+1,nxl:nxr,3) )
    769798    ALLOCATE ( dist_yz(nzb:nzt+1,nys:nyn,3) )
     
    831860    ELSE
    832861!
    833 !--    Define length scale for the imposed turbulence, which is defined as 8 times the minimum grid
    834 !--    spacing
    835        length_scale = 8.0_wp * MIN( dx, dy, MINVAL( dzw ) )
    836 !
    837 !--    Define constant to gradually decreasing length scales and Reynolds stress above ABL top.
    838 !--    Assure that no zero length scales are used.
    839        d_l = blend_coeff / MAX( length_scale, dx, dy, MINVAL( dzw ) )
     862!--    Precalulate the inversion number of xy-grid points - later used for mass conservation
     863       d_nxy = 1.0_wp / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
    840864!
    841865!--    Set flag indicating that turbulence is parametrized
     
    851875       CALL calc_scaling_variables
    852876!
    853 !--    Initialize lenght and time scales, which in turn are used to initialize the filter functions.
    854        CALL calc_length_and_time_scale
    855 !
    856877!--    Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'),
    857878!--    r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth,
    858879!--    friction velocity and velocity scale.
    859        CALL parametrize_reynolds_stress
     880       CALL parametrize_turbulence
    860881!
    861882!--    Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
     
    863884
    864885    ENDIF
    865 
    866 !
    867 !-- Assign initial profiles. Note, this is only required if turbulent inflow from the left is
    868 !-- desired, not in case of any of the nesting (offline or self nesting) approaches.
    869     IF ( .NOT. nesting_offline  .AND.  .NOT.  child_domain )  THEN
    870        u_init = mean_inflow_profiles(:,1)
    871        v_init = mean_inflow_profiles(:,2)
    872       !pt_init = mean_inflow_profiles(:,4)
    873        e_init = MAXVAL( mean_inflow_profiles(:,5) )
    874     ENDIF
    875 
    876 !
    877 !-- Define the size of the filter functions and allocate them.
    878     mergp = 0
    879 
    880     ! arrays must be large enough to cover the largest length scale
    881     DO  k = nzb, nzt+1
    882        j = MAX( ABS( nux(k) ), ABS( nuy(k) ), ABS( nuz(k) ),                                       &
    883                 ABS( nvx(k) ), ABS( nvy(k) ), ABS( nvz(k) ),                                       &
    884                 ABS( nwx(k) ), ABS( nwy(k) ), ABS( nwz(k) ) )
    885        IF ( j > mergp )  mergp = j
    886     ENDDO
    887 
    888 !     mergp  = 2 * mergp
    889 !     mergp = mergp
    890 
    891     ALLOCATE ( bux(-mergp:mergp,nzb:nzt+1),                                                        &
    892                buy(-mergp:mergp,nzb:nzt+1),                                                        &
    893                buz(-mergp:mergp,nzb:nzt+1),                                                        &
    894                bvx(-mergp:mergp,nzb:nzt+1),                                                        &
    895                bvy(-mergp:mergp,nzb:nzt+1),                                                        &
    896                bvz(-mergp:mergp,nzb:nzt+1),                                                        &
    897                bwx(-mergp:mergp,nzb:nzt+1),                                                        &
    898                bwy(-mergp:mergp,nzb:nzt+1),                                                        &
    899                bwz(-mergp:mergp,nzb:nzt+1) )
     886!
     887!-- Initial filter functions
     888    CALL stg_setup_filter_function
    900889
    901890!
     
    924913    fw_yz  = 0.0_wp
    925914    fwo_yz = 0.0_wp
    926 
    927 !
    928 !-- Create filter functions
    929     CALL stg_filter_func( nux, bux ) !filter ux
    930     CALL stg_filter_func( nuy, buy ) !filter uy
    931     CALL stg_filter_func( nuz, buz ) !filter uz
    932     CALL stg_filter_func( nvx, bvx ) !filter vx
    933     CALL stg_filter_func( nvy, bvy ) !filter vy
    934     CALL stg_filter_func( nvz, bvz ) !filter vz
    935     CALL stg_filter_func( nwx, bwx ) !filter wx
    936     CALL stg_filter_func( nwy, bwy ) !filter wy
    937     CALL stg_filter_func( nwz, bwz ) !filter wz
    938915
    939916#if defined( __parallel )
     
    10591036!-- core. In case there is only inflow from the left, it is sufficient to generate only random
    10601037!-- numbers for the yz-layer, else random numbers for the xz-layer are also required.
    1061     ALLOCATE ( id_rand_yz(-mergp+nys:nyn+mergp) )
    1062     ALLOCATE ( seq_rand_yz(5,-mergp+nys:nyn+mergp) )
     1038    ALLOCATE ( id_rand_yz(-mergp_limit+nys:nyn+mergp_limit) )
     1039    ALLOCATE ( seq_rand_yz(5,-mergp_limit+nys:nyn+mergp_limit) )
    10631040    id_rand_yz  = 0
    10641041    seq_rand_yz = 0
    10651042
    1066     CALL init_parallel_random_generator( ny, -mergp+nys, nyn+mergp, id_rand_yz, seq_rand_yz )
     1043    CALL init_parallel_random_generator( ny, -mergp_limit+nys, nyn+mergp_limit, id_rand_yz, seq_rand_yz )
    10671044
    10681045    IF ( nesting_offline  .OR.  ( child_domain .AND.  rans_mode_parent                             &
    10691046         .AND.  .NOT.  rans_mode ) )  THEN
    1070        ALLOCATE ( id_rand_xz(-mergp+nxl:nxr+mergp) )
    1071        ALLOCATE ( seq_rand_xz(5,-mergp+nxl:nxr+mergp) )
     1047       ALLOCATE ( id_rand_xz(-mergp_limit+nxl:nxr+mergp_limit) )
     1048       ALLOCATE ( seq_rand_xz(5,-mergp_limit+nxl:nxr+mergp_limit) )
    10721049       id_rand_xz  = 0
    10731050       seq_rand_xz = 0
    10741051
    1075        CALL init_parallel_random_generator( nx, -mergp+nxl, nxr+mergp, id_rand_xz, seq_rand_xz )
     1052       CALL init_parallel_random_generator( nx, -mergp_limit+nxl, nxr+mergp_limit, id_rand_xz, seq_rand_xz )
    10761053    ENDIF
    1077 
    1078 
    10791054
    10801055 END SUBROUTINE stg_init
     
    10861061!> Calculate filter function bxx from length scale nxx following Eg.9 and 10 (Xie and Castro, 2008)
    10871062!--------------------------------------------------------------------------------------------------!
    1088  SUBROUTINE stg_filter_func( nxx, bxx )
    1089 
    1090     INTEGER(iwp) ::  k    !< loop index
    1091     INTEGER(iwp) ::  n_k  !< length scale nXX in height k
    1092     INTEGER(iwp) ::  nf   !< index for length scales
     1063 SUBROUTINE stg_filter_func( nxx, bxx, mergp )
     1064
     1065    INTEGER(iwp) ::  k     !< loop index
     1066    INTEGER(iwp) ::  mergp !< passed length scale in grid points
     1067    INTEGER(iwp) ::  n_k   !< length scale nXX in height k
     1068    INTEGER(iwp) ::  nf    !< index for length scales
    10931069
    10941070    INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  nxx  !< length scale (in gp)
     
    11251101
    11261102
     1103!------------------------------------------------------------------------------!
     1104! Description:
     1105! ------------
     1106!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
     1107!> (Xie and Castro, 2008)
     1108!------------------------------------------------------------------------------!
     1109 SUBROUTINE stg_setup_filter_function
     1110
     1111    INTEGER(iwp) :: j  !< dummy value to calculate maximum length scale index
     1112    INTEGER(iwp) :: k  !< loop index along vertical direction
     1113!
     1114!-- Define the size of the filter functions and allocate them.
     1115    mergp_x = 0
     1116    mergp_y = 0
     1117    mergp_z = 0
     1118
     1119    DO  k = nzb, nzt+1
     1120       j = MAX( ABS(nux(k)), ABS(nvx(k)), ABS(nwx(k)) )
     1121       IF ( j > mergp_x )  mergp_x = j
     1122    ENDDO
     1123    DO  k = nzb, nzt+1
     1124       j = MAX( ABS(nuy(k)), ABS(nvy(k)), ABS(nwy(k)) )
     1125       IF ( j > mergp_y )  mergp_y = j
     1126    ENDDO
     1127    DO  k = nzb, nzt+1
     1128       j = MAX( ABS(nuz(k)), ABS(nvz(k)), ABS(nwz(k)) )
     1129       IF ( j > mergp_z )  mergp_z = j
     1130    ENDDO
     1131
     1132    mergp_xy = MAX( mergp_x, mergp_y )
     1133
     1134    IF ( ALLOCATED( bux ) )  DEALLOCATE( bux )
     1135    IF ( ALLOCATED( bvx ) )  DEALLOCATE( bvx )
     1136    IF ( ALLOCATED( bwx ) )  DEALLOCATE( bwx )
     1137    IF ( ALLOCATED( buy ) )  DEALLOCATE( buy )
     1138    IF ( ALLOCATED( bvy ) )  DEALLOCATE( bvy )
     1139    IF ( ALLOCATED( bwy ) )  DEALLOCATE( bwy )
     1140    IF ( ALLOCATED( buz ) )  DEALLOCATE( buz )
     1141    IF ( ALLOCATED( bvz ) )  DEALLOCATE( bvz )
     1142    IF ( ALLOCATED( bwz ) )  DEALLOCATE( bwz )
     1143
     1144    ALLOCATE ( bux(-mergp_x:mergp_x,nzb:nzt+1),                                                    &
     1145               buy(-mergp_y:mergp_y,nzb:nzt+1),                                                    &
     1146               buz(-mergp_z:mergp_z,nzb:nzt+1),                                                    &
     1147               bvx(-mergp_x:mergp_x,nzb:nzt+1),                                                    &
     1148               bvy(-mergp_y:mergp_y,nzb:nzt+1),                                                    &
     1149               bvz(-mergp_z:mergp_z,nzb:nzt+1),                                                    &
     1150               bwx(-mergp_x:mergp_x,nzb:nzt+1),                                                    &
     1151               bwy(-mergp_y:mergp_y,nzb:nzt+1),                                                    &
     1152               bwz(-mergp_z:mergp_z,nzb:nzt+1) )
     1153!
     1154!-- Create filter functions
     1155    CALL stg_filter_func( nux, bux, mergp_x ) !filter ux
     1156    CALL stg_filter_func( nuy, buy, mergp_y ) !filter uy
     1157    CALL stg_filter_func( nuz, buz, mergp_z ) !filter uz
     1158    CALL stg_filter_func( nvx, bvx, mergp_x ) !filter vx
     1159    CALL stg_filter_func( nvy, bvy, mergp_y ) !filter vy
     1160    CALL stg_filter_func( nvz, bvz, mergp_z ) !filter vz
     1161    CALL stg_filter_func( nwx, bwx, mergp_x ) !filter wx
     1162    CALL stg_filter_func( nwy, bwy, mergp_y ) !filter wy
     1163    CALL stg_filter_func( nwz, bwz, mergp_z ) !filter wz
     1164
     1165 END SUBROUTINE stg_setup_filter_function
     1166
    11271167!--------------------------------------------------------------------------------------------------!
    11281168! Description:
     
    12621302    LOGICAL  ::  stg_call = .FALSE.  !< control flag indicating whether turbulence was updated or only restored from last call
    12631303
    1264     REAL(wp) ::  dt_stg  !< wheighted subtimestep
     1304    REAL(wp) ::  dt_stg  !< time interval the STG is called
    12651305
    12661306    REAL(wp), DIMENSION(3) ::  mc_factor_l  !< local mass flux correction factor
     
    13441384!
    13451385!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
    1346                 IF ( tu(k) == 0.0_wp )  THEN
     1386                IF ( tu(k) == 0.0_wp  .OR.  adjustment_step )  THEN
    13471387                   fu_yz(k,j) = fuo_yz(k,j)
    13481388                ELSE
     
    13511391                ENDIF
    13521392
    1353                 IF ( tv(k) == 0.0_wp )  THEN
     1393                IF ( tv(k) == 0.0_wp  .OR.  adjustment_step )  THEN
    13541394                   fv_yz(k,j) = fvo_yz(k,j)
    13551395                ELSE
     
    13581398                ENDIF
    13591399
    1360                 IF ( tw(k) == 0.0_wp )  THEN
     1400                IF ( tw(k) == 0.0_wp  .OR.  adjustment_step )  THEN
    13611401                   fw_yz(k,j) = fwo_yz(k,j)
    13621402                ELSE
     
    13771417!
    13781418!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
    1379 !--             Additional factors are added to improve the variance of v and w
    1380                 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) *                              &
     1419                dist_yz(k,j,1) = a11(k) * fu_yz(k,j) *                                             &
    13811420                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
     1421
    13821422             ENDDO
    13831423          ENDDO
     
    13911431!--             Additional factors are added to improve the variance of v and w experimental test
    13921432!--             of 1.2
    1393                 dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL( a22 ) ) * 1.2_wp )                  &
    1394                                  * ( a21(k) * fu_yz(k,j) + a22(k) * fv_yz(k,j) ), 3.0_wp ) *       &
     1433                dist_yz(k,j,2) = ( a21(k) * fu_yz(k,j) + a22(k) * fv_yz(k,j) ) *                   &
    13951434                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    1396                 dist_yz(k,j,3) = MIN( ( SQRT( a33(k) / MAXVAL( a33 ) ) * 1.3_wp ) *                &
    1397                                       ( a31(k) * fu_yz(k,j) + a32(k) * fv_yz(k,j) + a33(k)         &
    1398                                         * fw_yz(k,j) ), 3.0_wp ) *                                 &
    1399                                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
     1435                dist_yz(k,j,3) = ( a31(k) * fu_yz(k,j) + a32(k) * fv_yz(k,j) +                     &
     1436                                   a33(k) * fw_yz(k,j) ) *                                         &
     1437                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
     1438
    14001439             ENDDO
    14011440          ENDDO
     
    15101549!
    15111550!--             Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
    1512                 IF ( tu(k) == 0.0_wp )  THEN
     1551                IF ( tu(k) == 0.0_wp  .OR.  adjustment_step )  THEN
    15131552                   fu_xz(k,i) = fuo_xz(k,i)
    15141553                ELSE
     
    15171556                ENDIF
    15181557
    1519                 IF ( tv(k) == 0.0_wp )  THEN
     1558                IF ( tv(k) == 0.0_wp  .OR.  adjustment_step )  THEN
    15201559                   fv_xz(k,i) = fvo_xz(k,i)
    15211560                ELSE
     
    15241563                ENDIF
    15251564
    1526                 IF ( tw(k) == 0.0_wp )  THEN
     1565                IF ( tw(k) == 0.0_wp  .OR.  adjustment_step )  THEN
    15271566                   fw_xz(k,i) = fwo_xz(k,i)
    15281567                ELSE
     
    15461585!--             Additional factors are added to improve the variance of v and w
    15471586                !experimental test of 1.2
    1548                 dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL( a22 ) ) * 1.2_wp )                  &
    1549                                  * ( a21(k) * fu_xz(k,i) + a22(k) * fv_xz(k,i) ), 3.0_wp ) *       &
     1587                dist_xz(k,i,2) = ( a21(k) * fu_xz(k,i) + a22(k) * fv_xz(k,i) ) *                   &
    15501588                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) )
    15511589             ENDDO
     
    15591597!--             Lund rotation following Eq. 17 in Xie and Castro (2008).
    15601598!--             Additional factors are added to improve the variance of v and w
    1561                 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) *                              &
     1599                dist_xz(k,i,1) = a11(k) * fu_xz(k,i) *                                             &
    15621600                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) )
    15631601
    1564                 dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL( a33 ) ) * 1.3_wp )                   &
    1565                                  * ( a31(k) * fu_xz(k,i) + a32(k) * fv_xz(k,i) + a33(k)            &
    1566                                  * fw_xz(k,i) ), 3.0_wp ) *                                        &
     1602                dist_xz(k,i,3) = ( a31(k) * fu_xz(k,i) + a32(k) * fv_xz(k,i) +                     &
     1603                                   a33(k) * fw_xz(k,i) ) *                                         &
    15671604                                 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) )
    15681605             ENDDO
     
    16541691!-- Finally, set time counter for calling STG to zero
    16551692    IF ( stg_call )  time_stg_call = 0.0_wp
     1693!
     1694!-- Set adjustment step to False to indicate that time correlation can be
     1695!-- switched-on again.
     1696    adjustment_step = .FALSE.
    16561697
    16571698    IF ( debug_output_timestep )  CALL debug_message( 'stg_main', 'end' )
     
    16841725    REAL(wp) ::  rand_sigma_inv  !< inverse of stdev of random number
    16851726
    1686     REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_y  !< filter function in y-direction
    1687     REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_z  !< filter function in z-direction
     1727    REAL(wp), DIMENSION(-mergp_y:mergp_y,nzb:nzt+1) :: b_y     !< filter function in y-direction
     1728    REAL(wp), DIMENSION(-mergp_z:mergp_z,nzb:nzt+1) :: b_z     !< filter function in z-direction
    16881729
    16891730    REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nys:nyn) ::  f_n_l  !<  local velocity seed
     
    16961737!-- must be defined globally in order to compute the correlation matrices. However, the calculation
    16971738!-- and normalization of random numbers is done locally, while the result is later distributed to
    1698 !-- all processes. Further, please note, a set of random numbers is only calculated for the left
    1699 !-- boundary, while the right boundary uses the same random numbers and thus also computes the same
     1739!-- all processes. Further, please note, a set of random numbers is only calculated for the west
     1740!-- boundary, while the east boundary uses the same random numbers and thus also computes the same
    17001741!-- correlation matrix.
    1701     ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nys:nyn+mergp) )
     1742    ALLOCATE( rand_it(nzb-mergp_z:nzt+1+mergp_z,-mergp_y+nys:nyn+mergp_y) )
    17021743    rand_it = 0.0_wp
    17031744
    17041745    rand_av        = 0.0_wp
    17051746    rand_sigma_inv = 0.0_wp
    1706     nyz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 )                      &
    1707                      * ( ny + mergp - ( 0 - mergp ) + 1 ), KIND = wp )
     1747    nyz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp_z - ( nzb - mergp_z ) + 1 )                  &
     1748                     * ( ny + mergp_y - ( 0 - mergp_y ) + 1 ), KIND = wp )
    17081749!
    17091750!-- Compute and normalize random numbers.
    1710     DO  j = nys - mergp, nyn + mergp
     1751    DO  j = nys - mergp_y, nyn + mergp_y
    17111752!
    17121753!--    Put the random seeds at grid point j
    17131754       CALL random_seed_parallel( put=seq_rand_yz(:,j) )
    1714        DO  k = nzb - mergp, nzt + 1 + mergp
     1755       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
    17151756          CALL random_number_parallel( random_dummy )
    17161757          rand_it(k,j) = random_dummy
     
    17241765!-- random numbers, the inner ghost layers mergp must not be summed-up, else the random numbers on
    17251766!-- the ghost layers will be stronger weighted as they also occur on the inner subdomains.
    1726     DO  j = MERGE( nys, nys - mergp, nys /= 0 ), MERGE( nyn, nyn + mergp, nyn /= ny )
    1727        DO  k = nzb - mergp, nzt + 1 + mergp
     1767    DO  j = MERGE( nys, nys - mergp_y, nys /= 0 ), MERGE( nyn, nyn + mergp_y, nyn /= ny )
     1768       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
    17281769          rand_av = rand_av + rand_it(k,j)
    17291770       ENDDO
    17301771    ENDDO
    1731 
     1772   
    17321773#if defined( __parallel )
    17331774!
     
    17411782!
    17421783!-- Now, compute the variance
    1743     DO  j = MERGE( nys, nys - mergp, nys /= 0 ), MERGE( nyn, nyn + mergp, nyn /= ny )
    1744        DO  k = nzb - mergp, nzt + 1 + mergp
     1784    DO  j = MERGE( nys, nys - mergp_y, nys /= 0 ), MERGE( nyn, nyn + mergp_y, nyn /= ny )
     1785       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
    17451786          rand_sigma_inv = rand_sigma_inv + rand_it(k,j)**2
    17461787       ENDDO
     
    18421883 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north )
    18431884
    1844     INTEGER(iwp) ::  i           !< loop index in x-direction
    1845     INTEGER(iwp) ::  id_north    !< core ids at respective boundaries
    1846     INTEGER(iwp) ::  id_south    !< core ids at respective boundaries
    1847     INTEGER(iwp) ::  ii          !< loop index in x-direction
    1848     INTEGER(iwp) ::  k           !< loop index in z-direction
    1849     INTEGER(iwp) ::  kk          !< loop index in z-direction
    1850     INTEGER(iwp) ::  send_count  !< send count for MPI_GATHERV
    1851 
    1852     INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  n_x  !< length scale in x-direction
    1853     INTEGER(iwp), DIMENSION(nzb:nzt+1) ::  n_z  !< length scale in z-direction
    1854 
    1855     REAL(wp) ::  nxz_inv         !< inverse of number of grid points in xz-slice
    1856     REAL(wp) ::  rand_av         !< average of random number
    1857     REAL(wp) ::  rand_sigma_inv  !< inverse of stdev of random number
    1858 
    1859     REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_x  !< filter function in x-direction
    1860     REAL(wp), DIMENSION(-mergp:mergp,nzb:nzt+1) ::  b_z  !< filter function in z-direction
    1861 
    1862     REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) ::  f_n_l  !<  local velocity seed
    1863     REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr)             ::  f_n    !<  velocity seed
    1864 
    1865     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it  !< global array of random numbers
     1885    INTEGER(iwp) :: i           !< loop index in x-direction
     1886    INTEGER(iwp) :: id_north    !< core ids at respective boundaries
     1887    INTEGER(iwp) :: id_south    !< core ids at respective boundaries
     1888    INTEGER(iwp) :: ii          !< loop index in x-direction
     1889    INTEGER(iwp) :: k           !< loop index in z-direction
     1890    INTEGER(iwp) :: kk          !< loop index in z-direction
     1891    INTEGER(iwp) :: send_count  !< send count for MPI_GATHERV
     1892
     1893    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x    !< length scale in x-direction
     1894    INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z    !< length scale in z-direction
     1895
     1896    REAL(wp) :: nxz_inv         !< inverse of number of grid points in xz-slice
     1897    REAL(wp) :: rand_av         !< average of random number
     1898    REAL(wp) :: rand_sigma_inv  !< inverse of stdev of random number
     1899
     1900    REAL(wp), DIMENSION(-mergp_x:mergp_x,nzb:nzt+1) :: b_x     !< filter function in x-direction
     1901    REAL(wp), DIMENSION(-mergp_z:mergp_z,nzb:nzt+1) :: b_z     !< filter function in z-direction
     1902   
     1903    REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) :: f_n_l   !<  local velocity seed
     1904    REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr)             :: f_n     !<  velocity seed
     1905
     1906    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  rand_it   !< global array of random numbers
    18661907
    18671908!
     
    18701911!-- must be defined globally in order to compute the correlation matrices. However, the calculation
    18711912!-- and normalization of random numbers is done locally, while the result is later distributed to
    1872 !-- all processes. Further, please note, a set of random numbers is only calculated for the left
    1873 !-- boundary, while the right boundary uses the same random numbers and thus also computes the same
     1913!-- all processes. Further, please note, a set of random numbers is only calculated for the south
     1914!-- boundary, while the north boundary uses the same random numbers and thus also computes the same
    18741915!-- correlation matrix.
    1875     ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp+nxl:nxr+mergp) )
     1916    ALLOCATE( rand_it(nzb-mergp_z:nzt+1+mergp_z,-mergp_x+nxl:nxr+mergp_x) )
    18761917    rand_it = 0.0_wp
    18771918
    18781919    rand_av        = 0.0_wp
    18791920    rand_sigma_inv = 0.0_wp
    1880     nxz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp - ( nzb - mergp ) + 1 )                      &
    1881                      * ( nx + mergp - ( 0 - mergp ) +1 ), KIND = wp )
     1921    nxz_inv        = 1.0_wp / REAL( ( nzt + 1 + mergp_z - ( nzb - mergp_z ) + 1 )                  &
     1922                     * ( nx + mergp_x - ( 0 - mergp_x ) +1 ), KIND = wp )
    18821923!
    18831924!-- Compute and normalize random numbers.
    1884     DO  i = nxl - mergp, nxr + mergp
     1925    DO  i = nxl - mergp_x, nxr + mergp_x
    18851926!
    18861927!--    Put the random seeds at grid point ii
    18871928       CALL random_seed_parallel( put=seq_rand_xz(:,i) )
    1888        DO  k = nzb - mergp, nzt + 1 + mergp
     1929       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
    18891930          CALL random_number_parallel( random_dummy )
    18901931          rand_it(k,i) = random_dummy
     
    18991940!-- summed-up, else the random numbers on the ghost layers will be stronger weighted as they
    19001941!-- also occur on the inner subdomains.
    1901     DO  i = MERGE( nxl, nxl - mergp, nxl /= 0 ), MERGE( nxr, nxr + mergp, nxr /= nx )
    1902        DO  k = nzb - mergp, nzt + 1 + mergp
     1942    DO  i = MERGE( nxl, nxl - mergp_x, nxl /= 0 ), MERGE( nxr, nxr + mergp_x, nxr /= nx )
     1943       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
    19031944          rand_av = rand_av + rand_it(k,i)
    19041945       ENDDO
    19051946    ENDDO
    1906 
     1947   
    19071948#if defined( __parallel )
    19081949!
     
    19161957!
    19171958!-- Now, compute the variance
    1918     DO  i = MERGE( nxl, nxl - mergp, nxl /= 0 ), MERGE( nxr, nxr + mergp, nxr /= nx )
    1919        DO  k = nzb - mergp, nzt + 1 + mergp
     1959    DO  i = MERGE( nxl, nxl - mergp_x, nxl /= 0 ), MERGE( nxr, nxr + mergp_x, nxr /= nx )
     1960       DO  k = nzb - mergp_z, nzt + 1 + mergp_z
    19201961          rand_sigma_inv = rand_sigma_inv + rand_it(k,i)**2
    19211962       ENDDO
     
    20042045! Description:
    20052046! ------------
    2006 !> Parametrization of the Reynolds stress tensor, following the parametrization described in Rotach
    2007 !> et al. (1996), which is applied in state-of-the-art dispserion modelling. Please note, the
    2008 !> parametrization does not distinguish between along-wind and cross-wind turbulence.
    2009 !--------------------------------------------------------------------------------------------------!
    2010  SUBROUTINE parametrize_reynolds_stress
    2011 
    2012     INTEGER(iwp) ::  k    !< loop index in z-direction
    2013 
    2014     REAL(wp)     ::  zzi  !< ratio of z/zi
    2015 
    2016     DO  k = nzb+1, nzt+1
     2047!> Parametrization of the Reynolds-stress componentes, turbulent length- and time scales. The
     2048!> parametrization follows Brost et al. (1982) with modifications described in Rotach et al. (1996),
     2049!> which is applied in state-of-the-art dispserion modelling.
     2050!--------------------------------------------------------------------------------------------------!
     2051  SUBROUTINE parametrize_turbulence
     2052
     2053    INTEGER(iwp) :: k                 !< loop index in z-direction
     2054
     2055    REAL(wp) ::  corr_term_uh         !< correction term in parametrization of horizontal variances for unstable stratification
     2056    REAL(wp) ::  d_zi                 !< inverse boundary-layer depth
     2057    REAL(wp) ::  length_scale_lon     !< length scale in flow direction
     2058    REAL(wp) ::  length_scale_lon_zi  !< length scale in flow direction at boundary-layer top
     2059    REAL(wp) ::  length_scale_lat     !< length scale in crosswind direction
     2060    REAL(wp) ::  length_scale_lat_zi  !< length scale in crosswind direction at boundary-layer top
     2061    REAL(wp) ::  length_scale_vert    !< length scale in vertical direction
     2062    REAL(wp) ::  length_scale_vert_zi !< length scale in vertical direction at boundary-layer top
     2063    REAL(wp) ::  time_scale_zi        !< time scale at boundary-layer top
     2064    REAL(wp) ::  r11_zi               !< longitudinal variance at boundary-layer top
     2065    REAL(wp) ::  r22_zi               !< crosswind variance at boundary-layer top
     2066    REAL(wp) ::  r33_zi               !< vertical variance at boundary-layer top
     2067    REAL(wp) ::  r31_zi               !< r31 at boundary-layer top
     2068    REAL(wp) ::  r32_zi               !< r32 at boundary-layer top
     2069    REAL(wp) ::  rlat1                !< first dummy argument for crosswind compoment of reynolds stress
     2070    REAL(wp) ::  rlat2                !< second dummy argument for crosswind compoment of reynolds stress
     2071    REAL(wp) ::  rlon1                !< first dummy argument for longitudinal compoment of reynolds stress
     2072    REAL(wp) ::  rlon2                !< second dummy argument for longitudinal compoment of reynolds stress
     2073    REAL(wp) ::  zzi                  !< ratio of z/zi
     2074
     2075    REAL(wp), DIMENSION(nzb+1:nzt+1) ::  cos_phi    !< cosine of angle between mean flow direction and x-axis
     2076    REAL(wp), DIMENSION(nzb+1:nzt+1) ::  phi        !< angle between mean flow direction and x-axis
     2077    REAL(wp), DIMENSION(nzb+1:nzt+1) ::  sin_phi    !< sine of angle between mean flow direction and x-axis
     2078
     2079!
     2080!-- Calculate the boundary-layer top index. Boundary-layer top is calculated by Richardson-bulk
     2081!-- criterion according to Heinze et al. (2017).
     2082    k_zi = MAX( 1, MINLOC( ABS( zu - zi_ribulk ), DIM = 1 ) )
     2083    IF ( zu(k_zi) > zi_ribulk )  k_zi = k_zi - 1
     2084    k_zi = MIN( k_zi, nzt )
     2085!
     2086!-- Calculate angle between flow direction and x-axis.
     2087    DO  k = nzb+1, nzt + 1
     2088       IF ( u_init(k) /= 0.0_wp )  THEN
     2089          phi(k) = ATAN( v_init(k) / u_init(k) )
     2090       ELSE
     2091          phi(k) = 0.5 * pi
     2092       ENDIF
     2093!
     2094!--    Pre-calculate sine and cosine.
     2095       cos_phi(k) = COS( phi(k) )
     2096       sin_phi(k) = SIN( phi(k) )
     2097    ENDDO
     2098!
     2099!-- Parametrize Reynolds-stress components. Please note, parametrization is formulated for the
     2100!-- stream- and spanwise components, while stream- and spanwise direction do not necessarily
     2101!-- coincide with the grid axis. Hence, components are rotated after computation.
     2102    d_zi = 1.0_wp / zi_ribulk
     2103    DO  k = nzb+1, k_zi
     2104
     2105       corr_term_uh = MERGE( 0.35_wp * ABS(                                                        &
     2106                        - zi_ribulk / ( kappa * scale_l - 10E-4_wp )                               &
     2107                                       )**( 2.0_wp / 3.0_wp ),                                     &
     2108                          0.0_wp,                                                                  &
     2109                          scale_l < 0.0_wp )
     2110!
     2111!--    Determine normalized height coordinate
     2112       zzi = zu(k) * d_zi
     2113!
     2114!--    Compute longitudinal and crosswind component of reynolds stress. Rotate these components by
     2115!--    the wind direction to map onto the u- and v-component. Note, since reynolds stress for
     2116!--    variances cannot become negative, take the absolute values.
     2117       rlon1 = scale_us**2 * ( corr_term_uh + 5.0_wp - 4.0_wp * zzi )
     2118       rlat1 = scale_us**2 * ( corr_term_uh + 2.0_wp - zzi )
     2119
     2120       r11(k) = ABS(   cos_phi(k) * rlon1 + sin_phi(k) * rlat1 )
     2121       r22(k) = ABS( - sin_phi(k) * rlon1 + cos_phi(k) * rlat1 )
     2122!
     2123!--    For u'v' no parametrization exist so far - ?. For simplicity assume a similar profile as
     2124!--    for u'w'.
     2125       r21(k) = 0.5_wp * ( r11(k) + r22(k) )
     2126!
     2127!--    w'w'
     2128       r33(k) = scale_wm**2 * (                                                                    &
     2129                   1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi )                        &
     2130                 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2                                   &
     2131                              )
     2132!
     2133!--    u'w' and v'w'. After calculation of the longitudinal and crosswind component
     2134!--    these are projected along the x- and y-direction. Note, it is assumed that
     2135!--    the flux within the boundary points opposite to the vertical gradient.
     2136       rlon2 = scale_us**2 * ( zzi - 1.0_wp )
     2137       rlat2 = scale_us**2 * ( 0.4 * zzi * ( 1.0_wp - zzi ) )
     2138
     2139       r31(k) = SIGN( ABS(   cos_phi(k) * rlon2 + sin_phi(k) * rlat2 ),                            &
     2140                   - ( u_init(k+1) - u_init(k-1) ) )
     2141       r32(k) = SIGN( ABS( - sin_phi(k) * rlon2 + cos_phi(k) * rlat2 ),                            &
     2142                   - ( v_init(k+1) - v_init(k-1) ) )
     2143!
     2144!--    Compute turbulent time scales according to Brost et al. (1982). Note, time scales are
     2145!--    limited to the adjustment time scales.
     2146       tu(k)  = MIN( dt_stg_adjust,                                                                &
     2147                     3.33_wp * zzi * ( 1.0 - 0.67_wp * zzi ) / scale_wm * zi_ribulk )
     2148
     2149       tv(k)  = tu(k)
     2150       tw(k)  = tu(k)
     2151
     2152       length_scale_lon  = MIN( SQRT( r11(k) ) * tu(k), zi_ribulk )
     2153       length_scale_lat  = MIN( SQRT( r22(k) ) * tv(k), zi_ribulk )
     2154       length_scale_vert = MIN( SQRT( r33(k) ) * tw(k), zi_ribulk )
     2155!
     2156!--    Assume isotropic turbulence length scales
     2157       nux(k) = MAX( INT( length_scale_lon  * ddx     ), 1 )
     2158       nuy(k) = MAX( INT( length_scale_lat  * ddy     ), 1 )
     2159       nvx(k) = MAX( INT( length_scale_lon  * ddx     ), 1 )
     2160       nvy(k) = MAX( INT( length_scale_lat  * ddy     ), 1 )
     2161       nwx(k) = MAX( INT( length_scale_lon  * ddx     ), 1 )
     2162       nwy(k) = MAX( INT( length_scale_lat  * ddy     ), 1 )
     2163       nuz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
     2164       nvz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
     2165       nwz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
     2166
     2167    ENDDO
     2168!
     2169!-- Above boundary-layer top length- and timescales as well as reynolds-stress components are
     2170!-- reduced to zero by a blending function.
     2171    length_scale_lon_zi  = SQRT( r11(k_zi) ) * tu(k_zi)
     2172    length_scale_lat_zi  = SQRT( r22(k_zi) ) * tv(k_zi)
     2173    length_scale_vert_zi = SQRT( r33(k_zi) ) * tu(k_zi)
     2174
     2175    time_scale_zi        = tu(k_zi)
     2176    r11_zi               = r11(k_zi)
     2177    r22_zi               = r22(k_zi)
     2178    r33_zi               = r33(k_zi)
     2179    r31_zi               = r31(k_zi)
     2180    r32_zi               = r32(k_zi)
     2181
     2182    d_l = blend_coeff / MAX( length_scale_vert_zi, dx, dy, MINVAL( dzw ) )
     2183
     2184    DO  k = k_zi+1, nzt+1
    20172185!
    20182186!--    Calculate function to gradually decrease Reynolds stress above ABL top.
    2019 !--    The function decreases to 1/10 after one length scale above the ABL top.
    20202187       blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) )
    2021 !
    2022 !--    Determine normalized height coordinate
    2023        zzi = zu(k) / zi_ribulk
    20242188!
    20252189!--    u'u' and v'v'. Assume isotropy. Note, add a small negative number to the denominator, else
    20262190!--    the mergpe-function can crash if scale_l is zero.
    2027        r11(k) = scale_us**2 * ( MERGE( 0.35_wp *                                                   &
    2028                          ABS( - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) )**( 2.0_wp / 3.0_wp ), &
    2029                          0.0_wp, scale_l < 0.0_wp ) + 5.0_wp - 4.0_wp * zzi                                 &
    2030                               ) * blend
    2031 
    2032        r22(k) = r11(k)
    2033 !
    2034 !--    w'w'
    2035        r33(k) = scale_wm**2 * ( 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi )           &
    2036                 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2 )  * blend
    2037 !
    2038 !--    u'w' and v'w'. Assume isotropy.
    2039        r31(k) = - scale_us**2 * ( 1.01_wp - MIN( zzi, 1.0_wp ) ) * blend
    2040        r32(k) = r31(k)
    2041 !
    2042 !--    For u'v' no parametrization exist so far - ?. For simplicity assume a similar profile as for
    2043 !--    u'w'.
    2044        r21(k) = r31(k)
     2191       r11(k) = r11_zi * blend
     2192       r22(k) = r22_zi * blend
     2193       r33(k) = r33_zi * blend
     2194       r31(k) = r31_zi * blend
     2195       r32(k) = r32_zi * blend
     2196       r21(k) = 0.5_wp * ( r11(k) + r22(k) )
     2197
     2198!
     2199!--    Compute turbulent time scales according to Brost et al. (1982).
     2200!--    Note, time scales are limited to the adjustment time scales.
     2201       tu(k)  = time_scale_zi * blend
     2202       tv(k)  = time_scale_zi * blend
     2203       tw(k)  = time_scale_zi * blend
     2204
     2205       length_scale_lon  = length_scale_lon_zi  * blend
     2206       length_scale_lat  = length_scale_lat_zi  * blend
     2207       length_scale_vert = length_scale_vert_zi * blend
     2208!
     2209!--    Assume isotropic turbulence length scales
     2210       nux(k) = MAX( INT( length_scale_lon  * ddx ),     1 )
     2211       nuy(k) = MAX( INT( length_scale_lat  * ddy ),     1 )
     2212       nvx(k) = MAX( INT( length_scale_lon  * ddx ),     1 )
     2213       nvy(k) = MAX( INT( length_scale_lat  * ddy ),     1 )
     2214       nwx(k) = MAX( INT( length_scale_lon  * ddx ),     1 )
     2215       nwy(k) = MAX( INT( length_scale_lat  * ddy ),     1 )
     2216       nuz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
     2217       nvz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
     2218       nwz(k) = MAX( INT( length_scale_vert * ddzw(k) ), 1 )
    20452219    ENDDO
    2046 
    2047 !
    2048 !-- Set bottom boundary condition
     2220!
     2221!-- Set bottom boundary condition for reynolds stresses
    20492222    r11(nzb) = r11(nzb+1)
    20502223    r22(nzb) = r22(nzb+1)
     
    20552228    r32(nzb) = r32(nzb+1)
    20562229
    2057 
    2058  END SUBROUTINE parametrize_reynolds_stress
     2230!
     2231!-- Set bottom boundary condition for the length and time scales
     2232    nux(nzb) = nux(nzb+1)
     2233    nuy(nzb) = nuy(nzb+1)
     2234    nuz(nzb) = nuz(nzb+1)
     2235    nvx(nzb) = nvx(nzb+1)
     2236    nvy(nzb) = nvy(nzb+1)
     2237    nvz(nzb) = nvz(nzb+1)
     2238    nwx(nzb) = nwx(nzb+1)
     2239    nwy(nzb) = nwy(nzb+1)
     2240    nwz(nzb) = nwz(nzb+1)
     2241
     2242    tu(nzb)  = tu(nzb+1)
     2243    tv(nzb)  = tv(nzb+1)
     2244    tw(nzb)  = tw(nzb+1)
     2245
     2246    adjustment_step = .TRUE.
     2247
     2248 END SUBROUTINE parametrize_turbulence
    20592249
    20602250!--------------------------------------------------------------------------------------------------!
     
    20712261    DO  k = nzb+1, nzt+1
    20722262       IF ( r11(k) > 10E-6_wp )  THEN
    2073           a11(k) = SQRT( r11(k) )
     2263          a11(k) = SQRT( r11(k) ) 
    20742264          a21(k) = r21(k) / a11(k)
    20752265          a31(k) = r31(k) / a11(k)
     
    20812271    ENDDO
    20822272    DO  k = nzb+1, nzt+1
    2083        a22(k) = r22(k) - a21(k)**2
     2273       a22(k) = r22(k) - a21(k)**2 
    20842274       IF ( a22(k) > 10E-6_wp )  THEN
    20852275          a22(k) = SQRT( a22(k) )
    2086           a32(k) = r32(k) - a21(k) * a31(k) / a22(k)
    2087        ELSE
     2276          a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
     2277       ELSE 
    20882278          a22(k) = 10E-8_wp
    20892279          a32(k) = 10E-8_wp
     
    21282318    CALL calc_scaling_variables
    21292319!
    2130 !-- Set length and time scales depending on boundary-layer height
    2131     CALL calc_length_and_time_scale
    2132 !
    2133 !-- Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'),
    2134 !-- r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth,
    2135 !-- friction velocity and velocity scale.
    2136     CALL parametrize_reynolds_stress
    2137 !
    2138 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
     2320!-- Parametrize Reynolds-stress tensor. Parametrization follows Brost et al. (1982) with
     2321!-- modifications described in Rotach et al. (1996) and is based on boundary-layer depth, friction
     2322!-- velocity and velocity scale.
     2323    CALL parametrize_turbulence
     2324!
     2325!-- Calculate coefficient matrix from Reynolds stress tensor 
     2326!-- (Lund rotation)
    21392327    CALL calc_coeff_matrix
    21402328!
    2141 !-- Determine filter functions on basis of updated length scales
    2142     CALL stg_filter_func( nux, bux ) !filter ux
    2143     CALL stg_filter_func( nuy, buy ) !filter uy
    2144     CALL stg_filter_func( nuz, buz ) !filter uz
    2145     CALL stg_filter_func( nvx, bvx ) !filter vx
    2146     CALL stg_filter_func( nvy, bvy ) !filter vy
    2147     CALL stg_filter_func( nvz, bvz ) !filter vz
    2148     CALL stg_filter_func( nwx, bwx ) !filter wx
    2149     CALL stg_filter_func( nwy, bwy ) !filter wy
    2150     CALL stg_filter_func( nwz, bwz ) !filter wz
     2329!-- Setup filter functions according to updated length scales.
     2330    CALL stg_setup_filter_function
    21512331!
    21522332!-- Reset time counter for controlling next adjustment to zero
     
    21572337 END SUBROUTINE stg_adjust
    21582338
    2159 
    2160 !--------------------------------------------------------------------------------------------------!
    2161 ! Description:
    2162 ! ------------
    2163 !> Calculates turbuluent length and time scales if these are not available from measurements.
    2164 !--------------------------------------------------------------------------------------------------!
    2165  SUBROUTINE calc_length_and_time_scale
    2166 
    2167     INTEGER(iwp) ::  k  !< loop index in z-direction
    2168 
    2169     REAL(wp) ::  length_scale_dum  !< effectively used length scale
    2170 
    2171 !
    2172 !-- In initial call the boundary-layer depth can be zero. This case, set minimum value for
    2173 !-- boundary-layer depth, to setup length scales correctly.
    2174     zi_ribulk = MAX( zi_ribulk, zw(nzb+2) )
    2175 !
    2176 !-- Set-up default turbulent length scales. From the numerical point of view the imposed
    2177 !-- perturbations should not be immediately dissipated by the numerics. The numerical dissipation,
    2178 !-- however, acts on scales up to 8 x the grid spacing. For this reason, set the turbulence length
    2179 !-- scale to 8 time the grid spacing. Further, above the boundary layer height, set turbulence
    2180 !-- lenght scales to zero (equivalent to not imposing any perturbations) in order to save
    2181 !-- computational costs. Typical time scales are derived by assuming Taylors's hypothesis, using
    2182 !-- the length scales and the mean profiles of u- and v-component.
    2183     DO  k = nzb+1, nzt+1
    2184 !
    2185 !--    Determine blending parameter. Within the boundary layer length scales are constant, while
    2186 !--    above lengths scales approach gradully zero, i.e. imposed turbulence is not correlated in
    2187 !--    space and time, just white noise, which saves computations power as the loops for the
    2188 !--    computation of the filter functions depend on the length scales. The value decreases to
    2189 !--    1/10 after one length scale above the ABL top.
    2190        blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) )
    2191 !
    2192 !--    Assume isotropic turbulence length scales
    2193        nux(k) = MAX( INT( length_scale * ddx ), 1 ) * blend
    2194        nuy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend
    2195        nvx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend
    2196        nvy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend
    2197        nwx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend
    2198        nwy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend
    2199 !
    2200 !--    Along the vertical direction limit the length scale further by the boundary-layer depth to
    2201 !--    assure that no length scales larger than the boundary-layer depth are used
    2202        length_scale_dum = MIN( length_scale, zi_ribulk )
    2203 
    2204        nuz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
    2205        nvz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
    2206        nwz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend
    2207 !
    2208 !--    Limit time scales, else they become very larger for low wind speed, imposing long-living
    2209 !--    inflow perturbations which in turn propagate further into the model domain. Use u_init and
    2210 !--    v_init to calculate the time scales, which will be equal to the inflow profiles, both, in
    2211 !--    offline nesting mode or in dirichlet/radiation mode.
    2212        tu(k) = MIN( dt_stg_adjust, length_scale / ( ABS( u_init(k) ) + 0.1_wp ) ) * blend
    2213        tv(k) = MIN( dt_stg_adjust, length_scale / ( ABS( v_init(k) ) + 0.1_wp ) ) * blend
    2214 !
    2215 !--    Time scale of w-component is a mixture from u- and v-component.
    2216        tw(k)  = SQRT( tu(k)**2 + tv(k)**2 ) * blend
    2217 
    2218     ENDDO
    2219 !
    2220 !-- Set bottom boundary condition for the length and time scales
    2221     nux(nzb) = nux(nzb+1)
    2222     nuy(nzb) = nuy(nzb+1)
    2223     nuz(nzb) = nuz(nzb+1)
    2224     nvx(nzb) = nvx(nzb+1)
    2225     nvy(nzb) = nvy(nzb+1)
    2226     nvz(nzb) = nvz(nzb+1)
    2227     nwx(nzb) = nwx(nzb+1)
    2228     nwy(nzb) = nwy(nzb+1)
    2229     nwz(nzb) = nwz(nzb+1)
    2230 
    2231     tu(nzb)  = tu(nzb+1)
    2232     tv(nzb)  = tv(nzb+1)
    2233     tw(nzb)  = tw(nzb+1)
    2234 
    2235 
    2236  END SUBROUTINE calc_length_and_time_scale
    22372339
    22382340!--------------------------------------------------------------------------------------------------!
     
    23052407#endif
    23062408
    2307     scale_us     = scale_us     / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
    2308     shf_mean     = shf_mean     / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
    2309     scale_l      = scale_l      / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
    2310     pt_surf_mean = pt_surf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp )
     2409    scale_us     = scale_us     * d_nxy
     2410    shf_mean     = shf_mean     * d_nxy
     2411    scale_l      = scale_l      * d_nxy
     2412    pt_surf_mean = pt_surf_mean * d_nxy
    23112413!
    23122414!-- Compute mean convective velocity scale. Note, in case the mean heat flux is negative, set
     
    23182420    ENDIF
    23192421!
     2422!-- At the initial run the friction velocity is initialized close to zero, leading to almost zero
     2423!-- disturbances at the boundaries. In order to obtain disturbances nevertheless, set a minium
     2424!-- value of friction velocity for the reynolds-stress parametrization.
     2425    IF ( time_since_reference_point <= 0.0_wp )  scale_us = MAX( scale_us, 0.05_wp )
     2426!
    23202427!-- Finally, in order to consider also neutral or stable stratification, compute momentum velocity
    23212428!-- scale from u* and convective velocity scale, according to Rotach et al. (1996).
Note: See TracChangeset for help on using the changeset viewer.