Changeset 4566
- Timestamp:
- Jun 16, 2020 10:11:51 AM (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r4562 r4566 25 25 ! ----------------- 26 26 ! $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 27 33 ! Parts of r4559 re-formatted 28 34 ! … … 151 157 !> parametrized within the boundary layer. 152 158 !> 153 !> @todo test restart 154 !> Enable cyclic_fill 159 !> @todo Enable cyclic_fill 155 160 !> Implement turbulence generation for e and pt 156 !> @todo Input of height-constant length scales via namelist157 161 !> @note <Enter notes on the module> 158 162 !> @bug Height information from input file is not used. Profiles from input must match with current … … 287 291 INTEGER(iwp) :: id_stg_right !< right lateral boundary core id in case of turbulence generator 288 292 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) 290 299 INTEGER(iwp) :: nzb_x_stg !< lower bound of z coordinate (required for transposing z on PEs along x) 291 300 INTEGER(iwp) :: nzt_x_stg !< upper bound of z coordinate (required for transposing z on PEs along x) … … 326 335 327 336 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 328 339 LOGICAL :: compute_velocity_seeds_local = .TRUE. !< switch to decide whether velocity seeds are computed locally 329 340 !< or if computation is distributed over several processes … … 336 347 REAL(wp) :: blend !< value to create gradually and smooth blending of Reynolds stress and length 337 348 !< scales above the boundary layer 338 REAL(wp) :: blend_coeff = - 2.3_wp !< coefficient used to ensure that blending functions decreases to 1/10 after349 REAL(wp) :: blend_coeff = -9.3_wp !< coefficient used to ensure that blending functions decreases to 1/10 after 339 350 !< one length scale above ABL top 340 351 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 342 354 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 spacing344 355 REAL(wp) :: scale_l !< scaling parameter used for turbulence parametrization - Obukhov length 345 356 REAL(wp) :: scale_us !< scaling parameter used for turbulence parametrization - friction velocity … … 390 401 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz !< velocity seed for w at yz plane with new random number 391 402 392 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_xz !< imposed disturbancesat north/south boundary393 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_yz !< imposed disturbancesat north/south boundary403 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 394 405 395 406 ! … … 756 767 ENDIF 757 768 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 758 777 759 778 ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), & … … 766 785 tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) ) 767 786 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 768 797 ALLOCATE ( dist_xz(nzb:nzt+1,nxl:nxr,3) ) 769 798 ALLOCATE ( dist_yz(nzb:nzt+1,nys:nyn,3) ) … … 831 860 ELSE 832 861 ! 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 ) 840 864 ! 841 865 !-- Set flag indicating that turbulence is parametrized … … 851 875 CALL calc_scaling_variables 852 876 ! 853 !-- Initialize lenght and time scales, which in turn are used to initialize the filter functions.854 CALL calc_length_and_time_scale855 !856 877 !-- Parametrize Reynolds-stress tensor, diagonal elements as well as r21 (v'u'), r31 (w'u'), 857 878 !-- r32 (w'v'). Parametrization follows Rotach et al. (1996) and is based on boundary-layer depth, 858 879 !-- friction velocity and velocity scale. 859 CALL parametrize_ reynolds_stress880 CALL parametrize_turbulence 860 881 ! 861 882 !-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation) … … 863 884 864 885 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 900 889 901 890 ! … … 924 913 fw_yz = 0.0_wp 925 914 fwo_yz = 0.0_wp 926 927 !928 !-- Create filter functions929 CALL stg_filter_func( nux, bux ) !filter ux930 CALL stg_filter_func( nuy, buy ) !filter uy931 CALL stg_filter_func( nuz, buz ) !filter uz932 CALL stg_filter_func( nvx, bvx ) !filter vx933 CALL stg_filter_func( nvy, bvy ) !filter vy934 CALL stg_filter_func( nvz, bvz ) !filter vz935 CALL stg_filter_func( nwx, bwx ) !filter wx936 CALL stg_filter_func( nwy, bwy ) !filter wy937 CALL stg_filter_func( nwz, bwz ) !filter wz938 915 939 916 #if defined( __parallel ) … … 1059 1036 !-- core. In case there is only inflow from the left, it is sufficient to generate only random 1060 1037 !-- 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) ) 1063 1040 id_rand_yz = 0 1064 1041 seq_rand_yz = 0 1065 1042 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 ) 1067 1044 1068 1045 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1069 1046 .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) ) 1072 1049 id_rand_xz = 0 1073 1050 seq_rand_xz = 0 1074 1051 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 ) 1076 1053 ENDIF 1077 1078 1079 1054 1080 1055 END SUBROUTINE stg_init … … 1086 1061 !> Calculate filter function bxx from length scale nxx following Eg.9 and 10 (Xie and Castro, 2008) 1087 1062 !--------------------------------------------------------------------------------------------------! 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 1093 1069 1094 1070 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: nxx !< length scale (in gp) … … 1125 1101 1126 1102 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 1127 1167 !--------------------------------------------------------------------------------------------------! 1128 1168 ! Description: … … 1262 1302 LOGICAL :: stg_call = .FALSE. !< control flag indicating whether turbulence was updated or only restored from last call 1263 1303 1264 REAL(wp) :: dt_stg !< wheighted subtimestep1304 REAL(wp) :: dt_stg !< time interval the STG is called 1265 1305 1266 1306 REAL(wp), DIMENSION(3) :: mc_factor_l !< local mass flux correction factor … … 1344 1384 ! 1345 1385 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1346 IF ( tu(k) == 0.0_wp ) THEN1386 IF ( tu(k) == 0.0_wp .OR. adjustment_step ) THEN 1347 1387 fu_yz(k,j) = fuo_yz(k,j) 1348 1388 ELSE … … 1351 1391 ENDIF 1352 1392 1353 IF ( tv(k) == 0.0_wp ) THEN1393 IF ( tv(k) == 0.0_wp .OR. adjustment_step ) THEN 1354 1394 fv_yz(k,j) = fvo_yz(k,j) 1355 1395 ELSE … … 1358 1398 ENDIF 1359 1399 1360 IF ( tw(k) == 0.0_wp ) THEN1400 IF ( tw(k) == 0.0_wp .OR. adjustment_step ) THEN 1361 1401 fw_yz(k,j) = fwo_yz(k,j) 1362 1402 ELSE … … 1377 1417 ! 1378 1418 !-- 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) * & 1381 1420 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1421 1382 1422 ENDDO 1383 1423 ENDDO … … 1391 1431 !-- Additional factors are added to improve the variance of v and w experimental test 1392 1432 !-- 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) ) * & 1395 1434 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 1400 1439 ENDDO 1401 1440 ENDDO … … 1510 1549 ! 1511 1550 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1512 IF ( tu(k) == 0.0_wp ) THEN1551 IF ( tu(k) == 0.0_wp .OR. adjustment_step ) THEN 1513 1552 fu_xz(k,i) = fuo_xz(k,i) 1514 1553 ELSE … … 1517 1556 ENDIF 1518 1557 1519 IF ( tv(k) == 0.0_wp ) THEN1558 IF ( tv(k) == 0.0_wp .OR. adjustment_step ) THEN 1520 1559 fv_xz(k,i) = fvo_xz(k,i) 1521 1560 ELSE … … 1524 1563 ENDIF 1525 1564 1526 IF ( tw(k) == 0.0_wp ) THEN1565 IF ( tw(k) == 0.0_wp .OR. adjustment_step ) THEN 1527 1566 fw_xz(k,i) = fwo_xz(k,i) 1528 1567 ELSE … … 1546 1585 !-- Additional factors are added to improve the variance of v and w 1547 1586 !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) ) * & 1550 1588 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 2 ) ) 1551 1589 ENDDO … … 1559 1597 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1560 1598 !-- 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) * & 1562 1600 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 1 ) ) 1563 1601 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) ) * & 1567 1604 MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 3 ) ) 1568 1605 ENDDO … … 1654 1691 !-- Finally, set time counter for calling STG to zero 1655 1692 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. 1656 1697 1657 1698 IF ( debug_output_timestep ) CALL debug_message( 'stg_main', 'end' ) … … 1684 1725 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1685 1726 1686 REAL(wp), DIMENSION(-mergp :mergp,nzb:nzt+1) :: b_y!< filter function in y-direction1687 REAL(wp), DIMENSION(-mergp :mergp,nzb:nzt+1) :: b_z!< filter function in z-direction1727 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 1688 1729 1689 1730 REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nys:nyn) :: f_n_l !< local velocity seed … … 1696 1737 !-- must be defined globally in order to compute the correlation matrices. However, the calculation 1697 1738 !-- 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 left1699 !-- boundary, while the right boundary uses the same random numbers and thus also computes the same1739 !-- 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 1700 1741 !-- 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) ) 1702 1743 rand_it = 0.0_wp 1703 1744 1704 1745 rand_av = 0.0_wp 1705 1746 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 ) 1708 1749 ! 1709 1750 !-- Compute and normalize random numbers. 1710 DO j = nys - mergp , nyn + mergp1751 DO j = nys - mergp_y, nyn + mergp_y 1711 1752 ! 1712 1753 !-- Put the random seeds at grid point j 1713 1754 CALL random_seed_parallel( put=seq_rand_yz(:,j) ) 1714 DO k = nzb - mergp , nzt + 1 + mergp1755 DO k = nzb - mergp_z, nzt + 1 + mergp_z 1715 1756 CALL random_number_parallel( random_dummy ) 1716 1757 rand_it(k,j) = random_dummy … … 1724 1765 !-- random numbers, the inner ghost layers mergp must not be summed-up, else the random numbers on 1725 1766 !-- 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 + mergp1767 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 1728 1769 rand_av = rand_av + rand_it(k,j) 1729 1770 ENDDO 1730 1771 ENDDO 1731 1772 1732 1773 #if defined( __parallel ) 1733 1774 ! … … 1741 1782 ! 1742 1783 !-- 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 + mergp1784 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 1745 1786 rand_sigma_inv = rand_sigma_inv + rand_it(k,j)**2 1746 1787 ENDDO … … 1842 1883 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id_south, id_north ) 1843 1884 1844 INTEGER(iwp) :: 1845 INTEGER(iwp) :: 1846 INTEGER(iwp) :: 1847 INTEGER(iwp) :: 1848 INTEGER(iwp) :: 1849 INTEGER(iwp) :: 1850 INTEGER(iwp) :: 1851 1852 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x!< length scale in x-direction1853 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z!< length scale in z-direction1854 1855 REAL(wp) :: 1856 REAL(wp) :: 1857 REAL(wp) :: 1858 1859 REAL(wp), DIMENSION(-mergp :mergp,nzb:nzt+1) :: b_x!< filter function in x-direction1860 REAL(wp), DIMENSION(-mergp :mergp,nzb:nzt+1) :: b_z!< filter function in z-direction1861 1862 REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxl:nxr) :: f_n_l!< local velocity seed1863 REAL(wp), DIMENSION(nzb:nzt+1,nxl:nxr) :: f_n!< velocity seed1864 1865 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< global array of random numbers1885 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 1866 1907 1867 1908 ! … … 1870 1911 !-- must be defined globally in order to compute the correlation matrices. However, the calculation 1871 1912 !-- 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 left1873 !-- boundary, while the rightboundary uses the same random numbers and thus also computes the same1913 !-- 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 1874 1915 !-- 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) ) 1876 1917 rand_it = 0.0_wp 1877 1918 1878 1919 rand_av = 0.0_wp 1879 1920 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 ) 1882 1923 ! 1883 1924 !-- Compute and normalize random numbers. 1884 DO i = nxl - mergp , nxr + mergp1925 DO i = nxl - mergp_x, nxr + mergp_x 1885 1926 ! 1886 1927 !-- Put the random seeds at grid point ii 1887 1928 CALL random_seed_parallel( put=seq_rand_xz(:,i) ) 1888 DO k = nzb - mergp , nzt + 1 + mergp1929 DO k = nzb - mergp_z, nzt + 1 + mergp_z 1889 1930 CALL random_number_parallel( random_dummy ) 1890 1931 rand_it(k,i) = random_dummy … … 1899 1940 !-- summed-up, else the random numbers on the ghost layers will be stronger weighted as they 1900 1941 !-- 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 + mergp1942 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 1903 1944 rand_av = rand_av + rand_it(k,i) 1904 1945 ENDDO 1905 1946 ENDDO 1906 1947 1907 1948 #if defined( __parallel ) 1908 1949 ! … … 1916 1957 ! 1917 1958 !-- 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 + mergp1959 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 1920 1961 rand_sigma_inv = rand_sigma_inv + rand_it(k,i)**2 1921 1962 ENDDO … … 2004 2045 ! Description: 2005 2046 ! ------------ 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 2017 2185 ! 2018 2186 !-- 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.2020 2187 blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) ) 2021 !2022 !-- Determine normalized height coordinate2023 zzi = zu(k) / zi_ribulk2024 2188 ! 2025 2189 !-- u'u' and v'v'. Assume isotropy. Note, add a small negative number to the denominator, else 2026 2190 !-- 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 ) 2045 2219 ENDDO 2046 2047 ! 2048 !-- Set bottom boundary condition 2220 ! 2221 !-- Set bottom boundary condition for reynolds stresses 2049 2222 r11(nzb) = r11(nzb+1) 2050 2223 r22(nzb) = r22(nzb+1) … … 2055 2228 r32(nzb) = r32(nzb+1) 2056 2229 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 2059 2249 2060 2250 !--------------------------------------------------------------------------------------------------! … … 2071 2261 DO k = nzb+1, nzt+1 2072 2262 IF ( r11(k) > 10E-6_wp ) THEN 2073 a11(k) = SQRT( r11(k) ) 2263 a11(k) = SQRT( r11(k) ) 2074 2264 a21(k) = r21(k) / a11(k) 2075 2265 a31(k) = r31(k) / a11(k) … … 2081 2271 ENDDO 2082 2272 DO k = nzb+1, nzt+1 2083 a22(k) = r22(k) - a21(k)**2 2273 a22(k) = r22(k) - a21(k)**2 2084 2274 IF ( a22(k) > 10E-6_wp ) THEN 2085 2275 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 2088 2278 a22(k) = 10E-8_wp 2089 2279 a32(k) = 10E-8_wp … … 2128 2318 CALL calc_scaling_variables 2129 2319 ! 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) 2139 2327 CALL calc_coeff_matrix 2140 2328 ! 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 2151 2331 ! 2152 2332 !-- Reset time counter for controlling next adjustment to zero … … 2157 2337 END SUBROUTINE stg_adjust 2158 2338 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_scale2166 2167 INTEGER(iwp) :: k !< loop index in z-direction2168 2169 REAL(wp) :: length_scale_dum !< effectively used length scale2170 2171 !2172 !-- In initial call the boundary-layer depth can be zero. This case, set minimum value for2173 !-- 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 imposed2177 !-- 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 length2179 !-- scale to 8 time the grid spacing. Further, above the boundary layer height, set turbulence2180 !-- lenght scales to zero (equivalent to not imposing any perturbations) in order to save2181 !-- computational costs. Typical time scales are derived by assuming Taylors's hypothesis, using2182 !-- the length scales and the mean profiles of u- and v-component.2183 DO k = nzb+1, nzt+12184 !2185 !-- Determine blending parameter. Within the boundary layer length scales are constant, while2186 !-- above lengths scales approach gradully zero, i.e. imposed turbulence is not correlated in2187 !-- space and time, just white noise, which saves computations power as the loops for the2188 !-- computation of the filter functions depend on the length scales. The value decreases to2189 !-- 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 scales2193 nux(k) = MAX( INT( length_scale * ddx ), 1 ) * blend2194 nuy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend2195 nvx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend2196 nvy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend2197 nwx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend2198 nwy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend2199 !2200 !-- Along the vertical direction limit the length scale further by the boundary-layer depth to2201 !-- assure that no length scales larger than the boundary-layer depth are used2202 length_scale_dum = MIN( length_scale, zi_ribulk )2203 2204 nuz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend2205 nvz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend2206 nwz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend2207 !2208 !-- Limit time scales, else they become very larger for low wind speed, imposing long-living2209 !-- inflow perturbations which in turn propagate further into the model domain. Use u_init and2210 !-- v_init to calculate the time scales, which will be equal to the inflow profiles, both, in2211 !-- offline nesting mode or in dirichlet/radiation mode.2212 tu(k) = MIN( dt_stg_adjust, length_scale / ( ABS( u_init(k) ) + 0.1_wp ) ) * blend2213 tv(k) = MIN( dt_stg_adjust, length_scale / ( ABS( v_init(k) ) + 0.1_wp ) ) * blend2214 !2215 !-- Time scale of w-component is a mixture from u- and v-component.2216 tw(k) = SQRT( tu(k)**2 + tv(k)**2 ) * blend2217 2218 ENDDO2219 !2220 !-- Set bottom boundary condition for the length and time scales2221 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_scale2237 2339 2238 2340 !--------------------------------------------------------------------------------------------------! … … 2305 2407 #endif 2306 2408 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 2311 2413 ! 2312 2414 !-- Compute mean convective velocity scale. Note, in case the mean heat flux is negative, set … … 2318 2420 ENDIF 2319 2421 ! 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 ! 2320 2427 !-- Finally, in order to consider also neutral or stable stratification, compute momentum velocity 2321 2428 !-- scale from u* and convective velocity scale, according to Rotach et al. (1996).
Note: See TracChangeset
for help on using the changeset viewer.
