Changeset 4022 for palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
- Timestamp:
- Jun 12, 2019 11:52:39 AM (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r3987 r4022 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Several bugfixes and improvements 28 ! - revise bias correction of the imposed perturbations (correction via volume 29 ! flow can create instabilities in case the mean volume flow is close to zero) 30 ! - introduce lower limits in calculation of coefficient matrix, else the 31 ! calculation may become numerically unstable 32 ! - impose perturbations every timestep, even though no new set of perturbations 33 ! is generated in case dt_stg_call /= dt_3d 34 ! - Implement a gradual decrease of Reynolds stress and length scales above 35 ! ABL height (within 1 length scale above ABL depth to 1/10) rather than a 36 ! discontinuous decrease 37 ! - Bugfix in non-nested case: use ABL height for parametrized turbulence 38 ! 39 ! 3987 2019-05-22 09:52:13Z kanani 27 40 ! Introduce alternative switch for debug output during timestepping 28 41 ! … … 206 219 207 220 USE arrays_3d, & 208 ONLY: mean_inflow_profiles, q, pt, u, v, w, zu, zw 221 ONLY: dzw, & 222 ddzw, & 223 drho_air, & 224 mean_inflow_profiles, & 225 q, & 226 pt, & 227 u, & 228 u_init, & 229 v, & 230 v_init, & 231 w, & 232 zu, & 233 zw 209 234 210 235 USE basic_constants_and_equations_mod, & 211 ONLY: g, kappa, pi 236 ONLY: g, & 237 kappa, & 238 pi 212 239 213 240 USE control_parameters, & 214 ONLY: debug_output_timestep, & 241 ONLY: bc_lr, & 242 bc_ns, & 243 child_domain, & 244 coupling_char, & 245 debug_output_timestep, & 246 dt_3d, & 247 e_init, & 215 248 initializing_actions, & 249 intermediate_timestep_count, & 250 intermediate_timestep_count_max, & 251 length, & 216 252 message_string, & 253 nesting_offline, & 217 254 num_mean_inflow_profiles, & 218 syn_turb_gen 255 rans_mode, & 256 restart_string, & 257 syn_turb_gen, & 258 time_since_reference_point, & 259 turbulent_inflow 260 261 USE grid_variables, & 262 ONLY: ddx, & 263 ddy, & 264 dx, & 265 dy 219 266 220 267 USE indices, & 221 ONLY: nbgp, nzb, nzt, nxl, nxlg, nxr, nxrg, nys, nyn, nyng, nysg 268 ONLY: nbgp, & 269 nz, & 270 nzb, & 271 nzt, & 272 nx, & 273 nxl, & 274 nxlu, & 275 nxlg, & 276 nxr, & 277 nxrg, & 278 ny, & 279 nys, & 280 nysv, & 281 nyn, & 282 nyng, & 283 nysg, & 284 wall_flags_0 222 285 223 286 USE kinds … … 228 291 229 292 USE nesting_offl_mod, & 230 ONLY: zi_ribulk 293 ONLY: nesting_offl_calc_zi, & 294 zi_ribulk 231 295 232 296 USE pegrid, & 233 ONLY: comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims 297 ONLY: comm1dx, & 298 comm1dy, & 299 comm2d, & 300 ierr, & 301 myidx, & 302 myidy, & 303 pdims 304 305 USE pmc_interface, & 306 ONLY : rans_mode_parent 234 307 235 308 USE transpose_indices, & 236 ONLY: nzb_x, nzt_x 309 ONLY: nzb_x, & 310 nzt_x 237 311 238 312 … … 263 337 INTEGER(iwp) :: nzt_y_stg !< upper bound of z coordinate (required for transposing z on PEs along y) 264 338 339 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_xz = 0 !< number of non-topography grid points at xz cross-sections, 340 !< required for bias correction of imposed perturbations 341 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_yz = 0 !< number of non-topography grid points at yz cross-sections, 342 !< required for bias correction of imposed perturbations 343 265 344 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz !< displacement for MPI_GATHERV 266 345 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz !< receive count for MPI_GATHERV … … 278 357 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed !< seed of random number for rn-generator 279 358 359 REAL(wp) :: blend !< value to create gradually and smooth blending of Reynolds stress and length 360 !< scales above the boundary layer 361 REAL(wp) :: blend_coeff = -2.3_wp !< coefficient used to ensure that blending functions decreases to 1/10 after 362 !< one length scale above ABL top 363 REAL(wp) :: d_l !< blend_coeff/length_scale 364 REAL(wp) :: length_scale !< length scale, default is 8 x minimum grid spacing 280 365 REAL(wp) :: dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics 281 366 REAL(wp) :: dt_stg_call = 5.0_wp !< time interval for calling synthetic turbulence generator 282 REAL(wp) :: mc_factor = 1.0_wp !< mass flux correction factor283 367 REAL(wp) :: scale_l !< scaling parameter used for turbulence parametrization - Obukhov length 284 368 REAL(wp) :: scale_us !< scaling parameter used for turbulence parametrization - friction velocity … … 286 370 REAL(wp) :: time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information 287 371 REAL(wp) :: time_stg_call = 0.0_wp !< time counter for calling generator 372 373 REAL(wp), DIMENSION(3) :: mc_factor = 1.0_wp !< correction factor for the u,v,w-components to maintain original mass flux 288 374 289 375 … … 421 507 SUBROUTINE stg_check_parameters 422 508 423 424 USE control_parameters, &425 ONLY: bc_lr, bc_ns, child_domain, nesting_offline, rans_mode, &426 turbulent_inflow427 428 USE pmc_interface, &429 ONLY : rans_mode_parent430 431 432 509 IMPLICIT NONE 433 510 … … 548 625 SUBROUTINE stg_init 549 626 550 551 USE arrays_3d, &552 ONLY: ddzw, u_init, v_init553 554 USE control_parameters, &555 ONLY: child_domain, coupling_char, e_init, nesting_offline, rans_mode556 557 USE grid_variables, &558 ONLY: ddy559 560 USE indices, &561 ONLY: nz562 563 USE pmc_interface, &564 ONLY : rans_mode_parent565 566 567 627 IMPLICIT NONE 568 628 … … 581 641 INTEGER(iwp) :: nseed !< dimension of random number seed 582 642 INTEGER(iwp) :: startseed = 1234567890 !< start seed for random number generator 583 643 644 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_xz_l = 0 !< number of non-topography grid points at xz-cross-section on subdomain 645 INTEGER(iwp), DIMENSION(3) :: nr_non_topo_yz_l = 0 !< number of non-topography grid points at yz-cross-section on subdomain 584 646 ! 585 647 !-- Dummy variables used for reading profiles … … 782 844 ELSE 783 845 ! 846 !-- Define length scale for the imposed turbulence, which is defined as 847 !-- 8 times the minimum grid spacing 848 length_scale = 8.0_wp * MIN( dx, dy, MINVAL( dzw ) ) 849 ! 850 !-- Define constant to gradually decrease length scales and Reynolds stress 851 !-- above ABL top. Assure that no zero length scales are used. 852 d_l = blend_coeff / MAX( length_scale, dx, dy, MINVAL( dzw ) ) 853 ! 784 854 !-- Set flag indicating that turbulence is parametrized 785 855 parametrize_inflow_turbulence = .TRUE. 856 ! 857 !-- In case of dirichlet inflow boundary conditions only at one lateral 858 !-- boundary, i.e. in the case no offline or self nesting is applied but 859 !-- synthetic turbulence shall be parametrized nevertheless, the 860 !-- boundary-layer depth need to determined first. 861 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) & 862 CALL nesting_offl_calc_zi 786 863 ! 787 864 !-- Determine boundary-layer depth, which is used to initialize lenght … … 897 974 DO k = nzb, nzt+1 898 975 IF ( a11(k) .NE. 0._wp ) THEN 899 fu_yz(k,j) = ( u(k,j,i) / mc_factor- u_init(k) ) / a11(k)976 fu_yz(k,j) = ( u(k,j,i) - u_init(k) ) / a11(k) 900 977 ELSE 901 978 fu_yz(k,j) = 0._wp … … 903 980 904 981 IF ( a22(k) .NE. 0._wp ) THEN 905 fv_yz(k,j) = ( v(k,j,i) / mc_factor - a21(k) * fu_yz(k,j) -&906 v_init(k) ) / a22(k)982 fv_yz(k,j) = ( v(k,j,i) - & 983 a21(k) * fu_yz(k,j) - v_init(k) ) / a22(k) 907 984 ELSE 908 985 fv_yz(k,j) = 0._wp … … 910 987 911 988 IF ( a33(k) .NE. 0._wp ) THEN 912 fw_yz(k,j) = ( w(k,j,i) / mc_factor - a31(k) * fu_yz(k,j) - & 913 a32(k) * fv_yz(k,j) ) / a33(k) 989 fw_yz(k,j) = ( w(k,j,i) - & 990 a31(k) * fu_yz(k,j) - a32(k) * & 991 fv_yz(k,j) ) / a33(k) 914 992 ELSE 915 993 fw_yz = 0._wp … … 929 1007 930 1008 IF ( a11(k) .NE. 0._wp ) THEN 931 fu_xz(k,i) = ( u(k,j,i) / mc_factor- u_init(k) ) / a11(k)1009 fu_xz(k,i) = ( u(k,j,i) - u_init(k) ) / a11(k) 932 1010 ELSE 933 1011 fu_xz(k,i) = 0._wp … … 935 1013 936 1014 IF ( a22(k) .NE. 0._wp ) THEN 937 fv_xz(k,i) = ( v(k,j,i) / mc_factor - a21(k) * fu_xz(k,i) -&938 v_init(k) ) / a22(k)1015 fv_xz(k,i) = ( v(k,j,i) - & 1016 a21(k) * fu_xz(k,i) - v_init(k) ) / a22(k) 939 1017 ELSE 940 1018 fv_xz(k,i) = 0._wp … … 942 1020 943 1021 IF ( a33(k) .NE. 0._wp ) THEN 944 fw_xz(k,i) = ( w(k,j,i) / mc_factor - a31(k) * fu_xz(k,i) - & 945 a32(k) * fv_xz(k,i) ) / a33(k) 1022 fw_xz(k,i) = ( w(k,j,i) - & 1023 a31(k) * fu_xz(k,i) - & 1024 a32(k) * fv_xz(k,i) ) / a33(k) 946 1025 ELSE 947 1026 fw_xz = 0._wp … … 952 1031 ENDIF 953 1032 ENDIF 1033 ! 1034 !-- Count the number of non-topography grid points at the boundaries where 1035 !-- perturbations are imposed. This number is later used for bias corrections 1036 !-- of the perturbations, i.e. to force that their mean is zero. Please note, 1037 !-- due to the asymetry of u and v along x and y direction, respectively, 1038 !-- different cases must be distinguished. 1039 IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN 1040 ! 1041 !-- Number of grid points where perturbations are imposed on u 1042 IF ( myidx == id_stg_left ) i = nxl 1043 IF ( myidx == id_stg_right ) i = nxr+1 1044 1045 nr_non_topo_yz_l(1) = SUM( & 1046 MERGE( 1, 0, & 1047 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 1 ) ) ) 1048 ! 1049 !-- Number of grid points where perturbations are imposed on v and w 1050 IF ( myidx == id_stg_left ) i = nxl-1 1051 IF ( myidx == id_stg_right ) i = nxr+1 1052 1053 nr_non_topo_yz_l(2) = SUM( & 1054 MERGE( 1, 0, & 1055 BTEST( wall_flags_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1056 nr_non_topo_yz_l(3) = SUM( & 1057 MERGE( 1, 0, & 1058 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1059 1060 #if defined( __parallel ) 1061 CALL MPI_ALLREDUCE( nr_non_topo_yz_l, nr_non_topo_yz, 3, MPI_INTEGER, & 1062 MPI_SUM, comm1dy, ierr ) 1063 #else 1064 nr_non_topo_yz = nr_non_topo_yz_l 1065 #endif 1066 ENDIF 1067 1068 IF ( myidy == id_stg_south .OR. myidy == id_stg_north ) THEN 1069 ! 1070 !-- Number of grid points where perturbations are imposed on v 1071 IF ( myidy == id_stg_south ) j = nys 1072 IF ( myidy == id_stg_north ) j = nyn+1 1073 1074 nr_non_topo_xz_l(2) = SUM( & 1075 MERGE( 1, 0, & 1076 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 2 ) ) ) 1077 ! 1078 !-- Number of grid points where perturbations are imposed on u and w 1079 IF ( myidy == id_stg_south ) j = nys-1 1080 IF ( myidy == id_stg_north ) j = nyn+1 1081 1082 nr_non_topo_xz_l(1) = SUM( & 1083 MERGE( 1, 0, & 1084 BTEST( wall_flags_0(nzb:nzt,j,nxlu:nxr), 1 ) ) ) 1085 nr_non_topo_xz_l(3) = SUM( & 1086 MERGE( 1, 0, & 1087 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 3 ) ) ) 1088 1089 #if defined( __parallel ) 1090 CALL MPI_ALLREDUCE( nr_non_topo_xz_l, nr_non_topo_xz, 3, MPI_INTEGER, & 1091 MPI_SUM, comm1dx, ierr ) 1092 #else 1093 nr_non_topo_xz = nr_non_topo_xz_l 1094 #endif 1095 ENDIF 1096 954 1097 955 1098 END SUBROUTINE stg_init … … 1059 1202 1060 1203 1061 USE control_parameters, &1062 ONLY: length, restart_string1063 1064 1065 1204 IMPLICIT NONE 1066 1205 … … 1072 1211 1073 1212 SELECT CASE ( restart_string(1:length) ) 1074 1075 CASE ( 'mc_factor' ) 1076 READ ( 13 ) mc_factor 1213 1214 CASE ( 'time_stg_adjust' ) 1215 READ ( 13 ) time_stg_adjust 1216 1217 CASE ( 'time_stg_call' ) 1218 READ ( 13 ) time_stg_call 1219 1077 1220 CASE ( 'use_syn_turb_gen' ) 1078 1221 READ ( 13 ) use_syn_turb_gen … … 1097 1240 1098 1241 IMPLICIT NONE 1099 1100 CALL wrd_write_string( 'mc_factor' ) 1101 WRITE ( 14 ) mc_factor 1242 1243 CALL wrd_write_string( 'time_stg_adjust' ) 1244 WRITE ( 14 ) time_stg_adjust 1245 1246 CALL wrd_write_string( 'time_stg_call' ) 1247 WRITE ( 14 ) time_stg_call 1102 1248 1103 1249 CALL wrd_write_string( 'use_syn_turb_gen' ) … … 1118 1264 SUBROUTINE stg_main 1119 1265 1120 1121 USE arrays_3d, &1122 ONLY: dzw1123 1124 USE control_parameters, &1125 ONLY: child_domain, dt_3d, &1126 nesting_offline, rans_mode, time_since_reference_point, &1127 volume_flow_initial1128 1129 USE grid_variables, &1130 ONLY: dx, dy1131 1132 USE indices, &1133 ONLY: wall_flags_01134 1135 USE pmc_interface, &1136 ONLY : rans_mode_parent1137 1138 1139 1266 IMPLICIT NONE 1140 1267 … … 1142 1269 INTEGER(iwp) :: j !< loop index in y-direction 1143 1270 INTEGER(iwp) :: k !< loop index in z-direction 1271 1272 LOGICAL :: stg_call = .FALSE. !< control flag indicating whether turbulence was updated or only restored from last call 1144 1273 1145 1274 REAL(wp) :: dt_stg !< wheighted subtimestep 1146 REAL(wp) :: mc_factor_l !< local mass flux correction factor 1147 REAL(wp) :: volume_flow !< mass flux through lateral boundary 1148 REAL(wp) :: volume_flow_l !< local mass flux through lateral boundary 1149 1275 1276 REAL(wp), DIMENSION(3) :: mc_factor_l !< local mass flux correction factor 1150 1277 1151 1278 IF ( debug_output_timestep ) CALL debug_message( 'stg_main', 'start' ) … … 1153 1280 !-- Calculate time step which is needed for filter functions 1154 1281 dt_stg = MAX( dt_3d, dt_stg_call ) 1282 ! 1283 !-- Check if synthetic turbulence generator needs to be called and new 1284 !-- perturbations need to be created or if old disturbances can be imposed 1285 !-- again. 1286 IF ( time_stg_call >= dt_stg_call .AND. & 1287 intermediate_timestep_count == intermediate_timestep_count_max ) THEN 1288 stg_call = .TRUE. 1289 ELSE 1290 stg_call = .FALSE. 1291 ENDIF 1155 1292 ! 1156 1293 !-- Initial value of fu, fv, fw … … 1183 1320 ! 1184 1321 !-- New set of fu, fv, fw 1185 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left ) 1186 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left ) 1187 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left ) 1188 1189 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1190 .AND. .NOT. rans_mode ) ) THEN 1191 ! 1192 !-- Generate turbulence at right boundary 1193 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right ) 1194 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right ) 1195 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right ) 1196 ! 1197 !-- Generate turbulence at north boundary 1198 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north ) 1199 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north ) 1200 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north ) 1201 ! 1202 !-- Generate turbulence at south boundary 1203 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south ) 1204 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south ) 1205 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south ) 1322 IF ( stg_call ) THEN 1323 ! 1324 !-- Generate turbulence at left boundary (required in all applications, 1325 !-- i.e. offline nesting and turbulent inflow from the left) 1326 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_left ) 1327 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left ) 1328 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left ) 1329 1330 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1331 .AND. .NOT. rans_mode ) ) THEN 1332 ! 1333 !-- Generate turbulence at right boundary 1334 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right ) 1335 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right ) 1336 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right ) 1337 ! 1338 !-- Generate turbulence at north boundary 1339 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north ) 1340 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north ) 1341 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north ) 1342 ! 1343 !-- Generate turbulence at south boundary 1344 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south ) 1345 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south ) 1346 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south ) 1347 ENDIF 1206 1348 ENDIF 1207 1349 … … 1209 1351 !-- Turbulence generation at left and or right boundary 1210 1352 IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN 1211 1212 DO j = nysg, nyng 1213 DO k = nzb, nzt + 1 1214 ! 1215 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1216 IF ( tu(k) == 0.0_wp ) THEN 1217 fu_yz(k,j) = fuo_yz(k,j) 1218 ELSE 1219 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + & 1220 fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) ) 1221 ENDIF 1222 1223 IF ( tv(k) == 0.0_wp ) THEN 1224 fv_yz(k,j) = fvo_yz(k,j) 1225 ELSE 1226 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + & 1227 fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) ) 1228 ENDIF 1229 1230 IF ( tw(k) == 0.0_wp ) THEN 1231 fw_yz(k,j) = fwo_yz(k,j) 1232 ELSE 1233 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + & 1234 fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) ) 1235 ENDIF 1236 ! 1237 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1238 !-- Additional factors are added to improve the variance of v and w 1239 IF( k == 0 ) THEN 1240 dist_yz(k,j,1) = 0.0_wp 1241 dist_yz(k,j,2) = 0.0_wp 1242 dist_yz(k,j,3) = 0.0_wp 1243 ELSE 1244 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) 1245 !experimental test of 1.2 1353 ! 1354 !-- Calculate new set of perturbations. Do this only at last RK3-substep and 1355 !-- when dt_stg_call is exceeded. Else the old set of perturbations is 1356 !-- imposed 1357 IF ( stg_call ) THEN 1358 1359 DO j = nysg, nyng 1360 DO k = nzb, nzt + 1 1361 ! 1362 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1363 IF ( tu(k) == 0.0_wp ) THEN 1364 fu_yz(k,j) = fuo_yz(k,j) 1365 ELSE 1366 fu_yz(k,j) = fu_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + & 1367 fuo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) ) 1368 ENDIF 1369 1370 IF ( tv(k) == 0.0_wp ) THEN 1371 fv_yz(k,j) = fvo_yz(k,j) 1372 ELSE 1373 fv_yz(k,j) = fv_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + & 1374 fvo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) ) 1375 ENDIF 1376 1377 IF ( tw(k) == 0.0_wp ) THEN 1378 fw_yz(k,j) = fwo_yz(k,j) 1379 ELSE 1380 fw_yz(k,j) = fw_yz(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + & 1381 fwo_yz(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) ) 1382 ENDIF 1383 ENDDO 1384 ENDDO 1385 1386 dist_yz(nzb,:,1) = 0.0_wp 1387 dist_yz(nzb,:,2) = 0.0_wp 1388 dist_yz(nzb,:,3) = 0.0_wp 1389 1390 IF ( myidx == id_stg_left ) i = nxl 1391 IF ( myidx == id_stg_right ) i = nxr+1 1392 DO j = nysg, nyng 1393 DO k = nzb+1, nzt + 1 1394 ! 1395 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1396 !-- Additional factors are added to improve the variance of v and w 1397 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) * & 1398 MERGE( 1.0_wp, 0.0_wp, & 1399 BTEST( wall_flags_0(k,j,i), 1 ) ) 1400 ENDDO 1401 ENDDO 1402 1403 IF ( myidx == id_stg_left ) i = nxl-1 1404 IF ( myidx == id_stg_right ) i = nxr+1 1405 DO j = nysg, nyng 1406 DO k = nzb+1, nzt + 1 1407 ! 1408 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1409 !-- Additional factors are added to improve the variance of v and w 1410 !-- experimental test of 1.2 1246 1411 dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) ) & 1247 1412 * 1.2_wp ) & 1248 1413 * ( a21(k) * fu_yz(k,j) & 1249 + a22(k) * fv_yz(k,j) ), 3.0_wp ) 1414 + a22(k) * fv_yz(k,j) ), 3.0_wp ) * & 1415 MERGE( 1.0_wp, 0.0_wp, & 1416 BTEST( wall_flags_0(k,j,i), 2 ) ) 1250 1417 dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) ) & 1251 1418 * 1.3_wp ) & 1252 1419 * ( a31(k) * fu_yz(k,j) & 1253 1420 + a32(k) * fv_yz(k,j) & 1254 + a33(k) * fw_yz(k,j) ), 3.0_wp ) 1255 ENDIF 1256 1421 + a33(k) * fw_yz(k,j) ), 3.0_wp ) * & 1422 MERGE( 1.0_wp, 0.0_wp, & 1423 BTEST( wall_flags_0(k,j,i), 3 ) ) 1424 ENDDO 1257 1425 ENDDO 1258 END DO1426 ENDIF 1259 1427 1260 1428 ! 1261 1429 !-- Mass flux correction following Kim et al. (2013) 1262 1430 !-- This correction factor insures that the mass flux is preserved at the 1263 !-- inflow boundary 1264 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 1265 mc_factor_l = 0.0_wp 1266 mc_factor = 0.0_wp 1267 DO j = nys, nyn 1268 DO k = nzb+1, nzt 1269 mc_factor_l = mc_factor_l + dzw(k) * & 1270 ( mean_inflow_profiles(k,1) + dist_yz(k,j,1) ) 1271 ENDDO 1272 ENDDO 1273 1274 #if defined( __parallel ) 1275 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM, & 1276 comm1dy, ierr ) 1277 #else 1278 mc_factor = mc_factor_l 1279 #endif 1280 1281 mc_factor = volume_flow_initial(1) / mc_factor 1282 1283 ! 1284 !-- Add disturbance at the inflow 1285 DO j = nysg, nyng 1286 DO k = nzb, nzt+1 1287 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + & 1288 dist_yz(k,j,1) ) * mc_factor 1289 v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + & 1290 dist_yz(k,j,2) ) * mc_factor 1291 w(k,j,-nbgp:-1) = dist_yz(k,j,3) * mc_factor 1292 ENDDO 1293 ENDDO 1294 1295 ELSE 1296 ! 1297 !-- First, calculate volume flow at yz boundary 1431 !-- inflow boundary. First, calculate mean value of the imposed 1432 !-- perturbations at yz boundary. 1433 !-- Note, this needs to be done only after the last RK3-substep, else 1434 !-- the boundary values won't be accessed. 1435 IF ( intermediate_timestep_count == intermediate_timestep_count_max ) THEN 1436 mc_factor_l = 0.0_wp 1437 mc_factor = 0.0_wp 1438 ! 1439 !-- Sum up the original volume flows (with and without perturbations). 1440 !-- Note, for non-normal components (here v and w) it is actually no 1441 !-- volume flow. 1298 1442 IF ( myidx == id_stg_left ) i = nxl 1299 1443 IF ( myidx == id_stg_right ) i = nxr+1 1300 1301 volume_flow_l = 0.0_wp 1302 volume_flow = 0.0_wp 1303 mc_factor_l = 0.0_wp 1304 mc_factor = 0.0_wp 1305 DO j = nys, nyn 1306 DO k = nzb+1, nzt 1307 volume_flow_l = volume_flow_l + u(k,j,i) * dzw(k) * dy & 1308 * MERGE( 1.0_wp, 0.0_wp, & 1309 BTEST( wall_flags_0(k,j,i), 1 ) ) 1310 1311 mc_factor_l = mc_factor_l + ( u(k,j,i) + dist_yz(k,j,1) ) & 1312 * dzw(k) * dy & 1313 * MERGE( 1.0_wp, 0.0_wp, & 1314 BTEST( wall_flags_0(k,j,i), 1 ) ) 1315 ENDDO 1316 ENDDO 1444 1445 mc_factor_l(1) = SUM( dist_yz(nzb:nzt,nys:nyn,1) & 1446 * MERGE( 1.0_wp, 0.0_wp, & 1447 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 1 ) ) ) 1448 1449 IF ( myidx == id_stg_left ) i = nxl-1 1450 IF ( myidx == id_stg_right ) i = nxr+1 1451 1452 mc_factor_l(2) = SUM( dist_yz(nzb:nzt,nysv:nyn,2) & 1453 * MERGE( 1.0_wp, 0.0_wp, & 1454 BTEST( wall_flags_0(nzb:nzt,nysv:nyn,i), 2 ) ) ) 1455 mc_factor_l(3) = SUM( dist_yz(nzb:nzt,nys:nyn,3) & 1456 * MERGE( 1.0_wp, 0.0_wp, & 1457 BTEST( wall_flags_0(nzb:nzt,nys:nyn,i), 3 ) ) ) 1458 1317 1459 #if defined( __parallel ) 1318 CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, &1319 1, MPI_REAL, MPI_SUM, comm1dy, ierr )1320 1460 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, & 1321 1, MPI_REAL, MPI_SUM, comm1dy, ierr ) 1322 #else 1323 volume_flow = volume_flow_l 1324 mc_factor = mc_factor_l 1461 3, MPI_REAL, MPI_SUM, comm1dy, ierr ) 1462 #else 1463 mc_factor = mc_factor_l 1325 1464 #endif 1326 1327 mc_factor = volume_flow / mc_factor 1328 1329 ! 1465 ! 1466 !-- Calculate correction factor and force zero mean perturbations. 1467 mc_factor = mc_factor / REAL( nr_non_topo_yz, KIND = wp ) 1468 1469 IF ( myidx == id_stg_left ) i = nxl 1470 IF ( myidx == id_stg_right ) i = nxr+1 1471 1472 dist_yz(:,:,1) = ( dist_yz(:,:,1) - mc_factor(1) ) & 1473 * MERGE( 1.0_wp, 0.0_wp, & 1474 BTEST( wall_flags_0(:,:,i), 1 ) ) 1475 1476 1477 IF ( myidx == id_stg_left ) i = nxl-1 1478 IF ( myidx == id_stg_right ) i = nxr+1 1479 1480 dist_yz(:,:,2) = ( dist_yz(:,:,2) - mc_factor(2) ) & 1481 * MERGE( 1.0_wp, 0.0_wp, & 1482 BTEST( wall_flags_0(:,:,i), 2 ) ) 1483 1484 dist_yz(:,:,3) = ( dist_yz(:,:,3) - mc_factor(3) ) & 1485 * MERGE( 1.0_wp, 0.0_wp, & 1486 BTEST( wall_flags_0(:,:,i), 3 ) ) 1487 ! 1330 1488 !-- Add disturbances 1331 1489 IF ( myidx == id_stg_left ) THEN 1332 DO j = nys, nyn 1333 DO k = nzb+1, nzt 1334 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) & 1335 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1336 BTEST( wall_flags_0(k,j,0), 1 ) ) 1337 u(k,j,-1) = u(k,j,0) 1338 v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) & 1339 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1340 BTEST( wall_flags_0(k,j,-1), 2 ) ) 1341 w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) & 1342 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1490 ! 1491 !-- For the left boundary distinguish between mesoscale offline / self 1492 !-- nesting and turbulent inflow at the left boundary only. In the latter 1493 !-- case turbulence is imposed on the mean inflow profiles. 1494 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 1495 ! 1496 !-- Add disturbance at the inflow 1497 DO j = nysg, nyng 1498 DO k = nzb, nzt+1 1499 u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + & 1500 dist_yz(k,j,1) ) & 1501 * MERGE( 1.0_wp, 0.0_wp, & 1502 BTEST( wall_flags_0(k,j,0), 1 ) ) 1503 v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + & 1504 dist_yz(k,j,2) ) & 1505 * MERGE( 1.0_wp, 0.0_wp, & 1506 BTEST( wall_flags_0(k,j,-1), 2 ) ) 1507 w(k,j,-nbgp:-1) = dist_yz(k,j,3) & 1508 * MERGE( 1.0_wp, 0.0_wp, & 1343 1509 BTEST( wall_flags_0(k,j,-1), 3 ) ) 1510 ENDDO 1344 1511 ENDDO 1345 ENDDO 1512 ELSE 1513 1514 DO j = nys, nyn 1515 DO k = nzb+1, nzt 1516 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) & 1517 * MERGE( 1.0_wp, 0.0_wp, & 1518 BTEST( wall_flags_0(k,j,0), 1 ) ) 1519 u(k,j,-1) = u(k,j,0) 1520 v(k,j,-1) = ( v(k,j,-1) + dist_yz(k,j,2) ) & 1521 * MERGE( 1.0_wp, 0.0_wp, & 1522 BTEST( wall_flags_0(k,j,-1), 2 ) ) 1523 w(k,j,-1) = ( w(k,j,-1) + dist_yz(k,j,3) ) & 1524 * MERGE( 1.0_wp, 0.0_wp, & 1525 BTEST( wall_flags_0(k,j,-1), 3 ) ) 1526 ENDDO 1527 ENDDO 1528 ENDIF 1346 1529 ENDIF 1347 1530 IF ( myidx == id_stg_right ) THEN … … 1349 1532 DO k = nzb+1, nzt 1350 1533 u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) & 1351 * mc_factor * MERGE( 1.0_wp, 0.0_wp,&1352 BTEST( wall_flags_0(k,j,nxr+1), 1 ) )1534 * MERGE( 1.0_wp, 0.0_wp, & 1535 BTEST( wall_flags_0(k,j,nxr+1), 1 ) ) 1353 1536 v(k,j,nxr+1) = ( v(k,j,nxr+1) + dist_yz(k,j,2) ) & 1354 * mc_factor * MERGE( 1.0_wp, 0.0_wp,&1355 BTEST( wall_flags_0(k,j,nxr+1), 2 ) )1537 * MERGE( 1.0_wp, 0.0_wp, & 1538 BTEST( wall_flags_0(k,j,nxr+1), 2 ) ) 1356 1539 w(k,j,nxr+1) = ( w(k,j,nxr+1) + dist_yz(k,j,3) ) & 1357 * mc_factor * MERGE( 1.0_wp, 0.0_wp,&1358 BTEST( wall_flags_0(k,j,nxr+1), 3 ) )1540 * MERGE( 1.0_wp, 0.0_wp, & 1541 BTEST( wall_flags_0(k,j,nxr+1), 3 ) ) 1359 1542 ENDDO 1360 1543 ENDDO 1361 1544 ENDIF 1362 1545 ENDIF 1363 1364 1546 ENDIF 1365 1547 ! 1366 1548 !-- Turbulence generation at north and south boundary 1367 1549 IF ( myidy == id_stg_north .OR. myidy == id_stg_south ) THEN 1368 1369 DO i = nxlg, nxrg 1370 DO k = nzb, nzt + 1 1371 ! 1372 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1373 IF ( tu(k) == 0.0_wp ) THEN 1374 fu_xz(k,i) = fuo_xz(k,i) 1375 ELSE 1376 fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + & 1377 fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) ) 1378 ENDIF 1379 1380 IF ( tv(k) == 0.0_wp ) THEN 1381 fv_xz(k,i) = fvo_xz(k,i) 1382 ELSE 1383 fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + & 1384 fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) ) 1385 ENDIF 1386 1387 IF ( tw(k) == 0.0_wp ) THEN 1388 fw_xz(k,i) = fwo_xz(k,i) 1389 ELSE 1390 fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + & 1391 fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) ) 1392 ENDIF 1393 ! 1394 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1395 !-- Additional factors are added to improve the variance of v and w 1396 IF( k == 0 ) THEN 1397 dist_xz(k,i,1) = 0.0_wp 1398 dist_xz(k,i,2) = 0.0_wp 1399 dist_xz(k,i,3) = 0.0_wp 1400 1401 ELSE 1402 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) 1403 !experimental test of 1.2 1550 ! 1551 !-- Calculate new set of perturbations. Do this only at last RK3-substep and 1552 !-- when dt_stg_call is exceeded. Else the old set of perturbations is 1553 !-- imposed 1554 IF ( stg_call ) THEN 1555 DO i = nxlg, nxrg 1556 DO k = nzb, nzt + 1 1557 ! 1558 !-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008) 1559 IF ( tu(k) == 0.0_wp ) THEN 1560 fu_xz(k,i) = fuo_xz(k,i) 1561 ELSE 1562 fu_xz(k,i) = fu_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + & 1563 fuo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) ) 1564 ENDIF 1565 1566 IF ( tv(k) == 0.0_wp ) THEN 1567 fv_xz(k,i) = fvo_xz(k,i) 1568 ELSE 1569 fv_xz(k,i) = fv_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + & 1570 fvo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) ) 1571 ENDIF 1572 1573 IF ( tw(k) == 0.0_wp ) THEN 1574 fw_xz(k,i) = fwo_xz(k,i) 1575 ELSE 1576 fw_xz(k,i) = fw_xz(k,i) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + & 1577 fwo_xz(k,i) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) ) 1578 ENDIF 1579 ENDDO 1580 ENDDO 1581 1582 1583 dist_xz(nzb,:,1) = 0.0_wp 1584 dist_xz(nzb,:,2) = 0.0_wp 1585 dist_xz(nzb,:,3) = 0.0_wp 1586 1587 IF ( myidy == id_stg_south ) j = nys 1588 IF ( myidy == id_stg_north ) j = nyn+1 1589 DO i = nxlg, nxrg 1590 DO k = nzb+1, nzt + 1 1591 ! 1592 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1593 !-- Additional factors are added to improve the variance of v and w 1594 !experimental test of 1.2 1404 1595 dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) ) & 1405 1596 * 1.2_wp ) & 1406 1597 * ( a21(k) * fu_xz(k,i) & 1407 + a22(k) * fv_xz(k,i) ), 3.0_wp ) 1598 + a22(k) * fv_xz(k,i) ), 3.0_wp ) * & 1599 MERGE( 1.0_wp, 0.0_wp, & 1600 BTEST( wall_flags_0(k,j,i), 2 ) ) 1601 ENDDO 1602 ENDDO 1603 1604 IF ( myidy == id_stg_south ) j = nys-1 1605 IF ( myidy == id_stg_north ) j = nyn+1 1606 DO i = nxlg, nxrg 1607 DO k = nzb+1, nzt + 1 1608 ! 1609 !-- Lund rotation following Eq. 17 in Xie and Castro (2008). 1610 !-- Additional factors are added to improve the variance of v and w 1611 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) * & 1612 MERGE( 1.0_wp, 0.0_wp, & 1613 BTEST( wall_flags_0(k,j,i), 1 ) ) 1614 1408 1615 dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) ) & 1409 1616 * 1.3_wp ) & 1410 1617 * ( a31(k) * fu_xz(k,i) & 1411 1618 + a32(k) * fv_xz(k,i) & 1412 + a33(k) * fw_xz(k,i) ), 3.0_wp ) 1413 ENDIF 1414 1415 ENDDO 1416 ENDDO 1417 1418 ! 1419 !-- Mass flux correction following Kim et al. (2013) 1420 !-- This correction factor insures that the mass flux is preserved at the 1421 !-- inflow boundary. 1422 !-- First, calculate volume flow at xz boundary 1423 IF ( myidy == id_stg_south ) j = nys 1424 IF ( myidy == id_stg_north ) j = nyn+1 1425 1426 volume_flow_l = 0.0_wp 1427 volume_flow = 0.0_wp 1428 mc_factor_l = 0.0_wp 1429 mc_factor = 0.0_wp 1430 DO i = nxl, nxr 1431 DO k = nzb+1, nzt 1432 volume_flow_l = volume_flow_l + v(k,j,i) * dzw(k) * dx & 1433 * MERGE( 1.0_wp, 0.0_wp, & 1434 BTEST( wall_flags_0(k,j,i), 2 ) ) 1435 1436 mc_factor_l = mc_factor_l + ( v(k,j,i) + dist_xz(k,i,2) ) & 1437 * dzw(k) * dx & 1438 * MERGE( 1.0_wp, 0.0_wp, & 1439 BTEST( wall_flags_0(k,j,i), 2 ) ) 1440 ENDDO 1441 ENDDO 1442 #if defined( __parallel ) 1443 CALL MPI_ALLREDUCE( volume_flow_l, volume_flow, & 1444 1, MPI_REAL, MPI_SUM, comm1dx, ierr ) 1445 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, & 1446 1, MPI_REAL, MPI_SUM, comm1dx, ierr ) 1447 #else 1448 volume_flow = volume_flow_l 1449 mc_factor = mc_factor_l 1450 #endif 1451 1452 mc_factor = volume_flow / mc_factor 1453 1454 ! 1455 !-- Add disturbances 1456 IF ( myidy == id_stg_south ) THEN 1457 1458 DO i = nxl, nxr 1459 DO k = nzb+1, nzt 1460 u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) ) & 1461 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1462 BTEST( wall_flags_0(k,-1,i), 1 ) ) 1463 v(k,0,i) = ( v(k,0,i) + dist_xz(k,i,2) ) & 1464 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1465 BTEST( wall_flags_0(k,0,i), 2 ) ) 1466 v(k,-1,i) = v(k,0,i) 1467 w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3) ) & 1468 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1469 BTEST( wall_flags_0(k,-1,i), 3 ) ) 1619 + a33(k) * fw_xz(k,i) ), 3.0_wp ) * & 1620 MERGE( 1.0_wp, 0.0_wp, & 1621 BTEST( wall_flags_0(k,j,i), 3 ) ) 1470 1622 ENDDO 1471 1623 ENDDO 1472 1624 ENDIF 1473 IF ( myidy == id_stg_north ) THEN 1474 1475 DO i = nxl, nxr 1476 DO k = nzb+1, nzt 1477 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) & 1478 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1479 BTEST( wall_flags_0(k,nyn+1,i), 1 ) ) 1480 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) & 1481 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1482 BTEST( wall_flags_0(k,nyn+1,i), 2 ) ) 1483 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) & 1484 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1485 BTEST( wall_flags_0(k,nyn+1,i), 3 ) ) 1625 1626 ! 1627 !-- Mass flux correction following Kim et al. (2013) 1628 !-- This correction factor insures that the mass flux is preserved at the 1629 !-- inflow boundary. First, calculate mean value of the imposed 1630 !-- perturbations at yz boundary. 1631 !-- Note, this needs to be done only after the last RK3-substep, else 1632 !-- the boundary values won't be accessed. 1633 IF ( intermediate_timestep_count == intermediate_timestep_count_max ) THEN 1634 mc_factor_l = 0.0_wp 1635 mc_factor = 0.0_wp 1636 1637 IF ( myidy == id_stg_south ) j = nys 1638 IF ( myidy == id_stg_north ) j = nyn+1 1639 1640 mc_factor_l(2) = SUM( dist_xz(nzb:nzt,nxl:nxr,2) & 1641 * MERGE( 1.0_wp, 0.0_wp, & 1642 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 2 ) ) ) 1643 1644 IF ( myidy == id_stg_south ) j = nys-1 1645 IF ( myidy == id_stg_north ) j = nyn+1 1646 1647 mc_factor_l(1) = SUM( dist_xz(nzb:nzt,nxlu:nxr,1) & 1648 * MERGE( 1.0_wp, 0.0_wp, & 1649 BTEST( wall_flags_0(nzb:nzt,j,nxlu:nxr), 1 ) ) ) 1650 mc_factor_l(3) = SUM( dist_xz(nzb:nzt,nxl:nxr,3) & 1651 * MERGE( 1.0_wp, 0.0_wp, & 1652 BTEST( wall_flags_0(nzb:nzt,j,nxl:nxr), 3 ) ) ) 1653 1654 #if defined( __parallel ) 1655 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, & 1656 3, MPI_REAL, MPI_SUM, comm1dx, ierr ) 1657 #else 1658 mc_factor = mc_factor_l 1659 #endif 1660 1661 mc_factor = mc_factor / REAL( nr_non_topo_xz, KIND = wp ) 1662 1663 IF ( myidy == id_stg_south ) j = nys 1664 IF ( myidy == id_stg_north ) j = nyn+1 1665 1666 dist_xz(:,:,2) = ( dist_xz(:,:,2) - mc_factor(2) ) & 1667 * MERGE( 1.0_wp, 0.0_wp, & 1668 BTEST( wall_flags_0(:,j,:), 2 ) ) 1669 1670 1671 IF ( myidy == id_stg_south ) j = nys-1 1672 IF ( myidy == id_stg_north ) j = nyn+1 1673 1674 dist_xz(:,:,1) = ( dist_xz(:,:,1) - mc_factor(1) ) & 1675 * MERGE( 1.0_wp, 0.0_wp, & 1676 BTEST( wall_flags_0(:,j,:), 1 ) ) 1677 1678 dist_xz(:,:,3) = ( dist_xz(:,:,3) - mc_factor(3) ) & 1679 * MERGE( 1.0_wp, 0.0_wp, & 1680 BTEST( wall_flags_0(:,j,:), 3 ) ) 1681 ! 1682 !-- Add disturbances 1683 IF ( myidy == id_stg_south ) THEN 1684 DO i = nxl, nxr 1685 DO k = nzb+1, nzt 1686 u(k,-1,i) = ( u(k,-1,i) + dist_xz(k,i,1) ) & 1687 * MERGE( 1.0_wp, 0.0_wp, & 1688 BTEST( wall_flags_0(k,-1,i), 1 ) ) 1689 v(k,0,i) = ( v(k,0,i) + dist_xz(k,i,2) ) & 1690 * MERGE( 1.0_wp, 0.0_wp, & 1691 BTEST( wall_flags_0(k,0,i), 2 ) ) 1692 v(k,-1,i) = v(k,0,i) 1693 w(k,-1,i) = ( w(k,-1,i) + dist_xz(k,i,3) ) & 1694 * MERGE( 1.0_wp, 0.0_wp, & 1695 BTEST( wall_flags_0(k,-1,i), 3 ) ) 1696 ENDDO 1486 1697 ENDDO 1487 ENDDO 1698 ENDIF 1699 IF ( myidy == id_stg_north ) THEN 1700 1701 DO i = nxl, nxr 1702 DO k = nzb+1, nzt 1703 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) & 1704 * MERGE( 1.0_wp, 0.0_wp, & 1705 BTEST( wall_flags_0(k,nyn+1,i), 1 ) ) 1706 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) & 1707 * MERGE( 1.0_wp, 0.0_wp, & 1708 BTEST( wall_flags_0(k,nyn+1,i), 2 ) ) 1709 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) & 1710 * MERGE( 1.0_wp, 0.0_wp, & 1711 BTEST( wall_flags_0(k,nyn+1,i), 3 ) ) 1712 ENDDO 1713 ENDDO 1714 ENDIF 1488 1715 ENDIF 1489 1716 ENDIF 1490 1717 ! 1491 1718 !-- Finally, set time counter for calling STG to zero 1492 time_stg_call = 0.0_wp1719 IF ( stg_call ) time_stg_call = 0.0_wp 1493 1720 1494 1721 IF ( debug_output_timestep ) CALL debug_message( 'stg_main', 'end' ) … … 1506 1733 !------------------------------------------------------------------------------! 1507 1734 SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n, id ) 1508 1509 1510 USE indices, &1511 ONLY: ny1512 1735 1513 1736 IMPLICIT NONE … … 1633 1856 SUBROUTINE stg_generate_seed_xz( n_x, n_z, b_x, b_z, f_n, id ) 1634 1857 1635 1636 USE indices, &1637 ONLY: nx1638 1639 1640 1858 IMPLICIT NONE 1641 1859 … … 1758 1976 SUBROUTINE parametrize_reynolds_stress 1759 1977 1760 USE arrays_3d, &1761 ONLY: zu1762 1763 1978 IMPLICIT NONE 1764 1979 … … 1767 1982 REAL(wp) :: zzi !< ratio of z/zi 1768 1983 1769 !1770 !--1771 1984 DO k = nzb+1, nzt+1 1772 1773 IF ( zu(k) <= zi_ribulk ) THEN 1774 ! 1775 !-- Determine normalized height coordinate 1776 zzi = zu(k) / zi_ribulk 1777 ! 1778 !-- u'u' and v'v'. Assume isotropy. Note, add a small negative number 1779 !-- to the denominator, else the merge-function can crash if scale_l is 1780 !-- zero. 1781 r11(k) = scale_us**2 * ( & 1782 MERGE( 0.35_wp * ABS( & 1783 - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) & 1784 )**( 2.0_wp / 3.0_wp ), & 1785 0.0_wp, & 1786 scale_l < 0.0_wp ) & 1787 + 5.0_wp - 4.0_wp * zzi & 1788 ) 1789 1790 r22(k) = r11(k) 1791 ! 1792 !-- w'w' 1793 r33(k) = scale_wm**2 * ( & 1794 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) & 1795 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2 & 1796 ) 1797 ! 1798 !-- u'w' and v'w'. Assume isotropy. 1799 r31(k) = - scale_us**2 * ( & 1800 1.0_wp - EXP( 3.0_wp * ( zzi - 1.0_wp ) ) & 1801 ) 1802 1803 r32(k) = r31(k) 1804 ! 1805 !-- For u'v' no parametrization exist so far - ?. For simplicity assume 1806 !-- a similar profile as for u'w'. 1807 r21(k) = r31(k) 1808 ! 1809 !-- Above the boundary layer, assmume laminar flow conditions. 1810 ELSE 1811 r11(k) = 10E-8_wp 1812 r22(k) = 10E-8_wp 1813 r33(k) = 10E-8_wp 1814 r21(k) = 10E-8_wp 1815 r31(k) = 10E-8_wp 1816 r32(k) = 10E-8_wp 1817 ENDIF 1985 ! 1986 !-- Calculate function to gradually decrease Reynolds stress above ABL top. 1987 !-- The function decreases to 1/10 after one length scale above the 1988 !-- ABL top. 1989 blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) ) 1990 ! 1991 !-- Determine normalized height coordinate 1992 zzi = zu(k) / zi_ribulk 1993 ! 1994 !-- u'u' and v'v'. Assume isotropy. Note, add a small negative number 1995 !-- to the denominator, else the merge-function can crash if scale_l is 1996 !-- zero. 1997 r11(k) = scale_us**2 * ( & 1998 MERGE( 0.35_wp * ABS( & 1999 - zi_ribulk / ( kappa * scale_l - 10E-4_wp ) & 2000 )**( 2.0_wp / 3.0_wp ), & 2001 0.0_wp, & 2002 scale_l < 0.0_wp ) & 2003 + 5.0_wp - 4.0_wp * zzi & 2004 ) * blend 2005 2006 r22(k) = r11(k) 2007 ! 2008 !-- w'w' 2009 r33(k) = scale_wm**2 * ( & 2010 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( -2.0_wp * zzi ) & 2011 + ( 1.7_wp - zzi ) * ( scale_us / scale_wm )**2 & 2012 ) * blend 2013 ! 2014 !-- u'w' and v'w'. Assume isotropy. 2015 r31(k) = - scale_us**2 * ( 1.01_wp - MIN( zzi, 1.0_wp ) ) * blend 2016 r32(k) = r31(k) 2017 ! 2018 !-- For u'v' no parametrization exist so far - ?. For simplicity assume 2019 !-- a similar profile as for u'w'. 2020 r21(k) = r31(k) 1818 2021 ENDDO 1819 2022 … … 1845 2048 !-- Calculate coefficient matrix. Split loops to allow for loop vectorization. 1846 2049 DO k = nzb+1, nzt+1 1847 IF ( r11(k) > 0.0_wp ) THEN2050 IF ( r11(k) > 10E-6_wp ) THEN 1848 2051 a11(k) = SQRT( r11(k) ) 1849 2052 a21(k) = r21(k) / a11(k) … … 1857 2060 DO k = nzb+1, nzt+1 1858 2061 a22(k) = r22(k) - a21(k)**2 1859 IF ( a22(k) > 0.0_wp ) THEN2062 IF ( a22(k) > 10E-6_wp ) THEN 1860 2063 a22(k) = SQRT( a22(k) ) 1861 2064 a32(k) = r32(k) - a21(k) * a31(k) / a22(k) … … 1867 2070 DO k = nzb+1, nzt+1 1868 2071 a33(k) = r33(k) - a31(k)**2 - a32(k)**2 1869 IF ( a33(k) > 0.0_wp ) THEN2072 IF ( a33(k) > 10E-6_wp ) THEN 1870 2073 a33(k) = SQRT( a33(k) ) 1871 2074 ELSE … … 1897 2100 IF ( debug_output_timestep ) CALL debug_message( 'stg_adjust', 'start' ) 1898 2101 ! 1899 !-- Compute mean boundary layer height according to Richardson-Bulk 1900 !-- criterion using the inflow profiles. Further velocity scale as well as 1901 !-- mean friction velocity are calculated. 2102 !-- In case of dirichlet inflow boundary conditions only at one lateral 2103 !-- boundary, i.e. in the case no offline or self nesting is applied but 2104 !-- synthetic turbulence shall be parametrized nevertheless, the 2105 !-- boundary-layer depth need to determined first. 2106 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) & 2107 CALL nesting_offl_calc_zi 2108 ! 2109 !-- Compute scaling parameters (domain-averaged), such as friction velocity 2110 !-- are calculated. 1902 2111 CALL calc_scaling_variables 1903 2112 ! … … 1942 2151 SUBROUTINE calc_length_and_time_scale 1943 2152 1944 USE arrays_3d, &1945 ONLY: dzw, ddzw, u_init, v_init1946 1947 USE grid_variables, &1948 ONLY: ddx, ddy, dx, dy1949 1950 2153 IMPLICIT NONE 1951 2154 1952 2155 1953 INTEGER(iwp) :: k!< loop index in z-direction1954 1955 REAL(wp) :: length_scale !< typicallength scale2156 INTEGER(iwp) :: k !< loop index in z-direction 2157 2158 REAL(wp) :: length_scale_dum !< effectively used length scale 1956 2159 1957 2160 ! … … 1970 2173 !-- using the length scales and the mean profiles of u- and v-component. 1971 2174 DO k = nzb+1, nzt+1 1972 1973 length_scale = 8.0_wp * MIN( dx, dy, dzw(k) ) 1974 1975 IF ( zu(k) <= zi_ribulk ) THEN 1976 ! 1977 !-- Assume isotropic turbulence length scales 1978 nux(k) = MAX( INT( length_scale * ddx ), 1 ) 1979 nuy(k) = MAX( INT( length_scale * ddy ), 1 ) 1980 nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 ) 1981 nvx(k) = MAX( INT( length_scale * ddx ), 1 ) 1982 nvy(k) = MAX( INT( length_scale * ddy ), 1 ) 1983 nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 ) 1984 nwx(k) = MAX( INT( length_scale * ddx ), 1 ) 1985 nwy(k) = MAX( INT( length_scale * ddy ), 1 ) 1986 nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 ) 1987 ! 1988 !-- Limit time scales, else they become very larger for low wind speed, 1989 !-- imposing long-living inflow perturbations which in turn propagate 1990 !-- further into the model domain. Use u_init and v_init to calculate 1991 !-- the time scales, which will be equal to the inflow profiles, both, 1992 !-- in offline nesting mode or in dirichlet/radiation mode. 1993 tu(k) = MIN( dt_stg_adjust, length_scale / & 1994 ( ABS( u_init(k) ) + 0.1_wp ) ) 1995 tv(k) = MIN( dt_stg_adjust, length_scale / & 1996 ( ABS( v_init(k) ) + 0.1_wp ) ) 1997 ! 1998 !-- Time scale of w-component is a mixture from u- and v-component. 1999 tw(k) = SQRT( tu(k)**2 + tv(k)**2 ) 2000 ! 2001 !-- Above the boundary layer length scales are zero, i.e. imposed turbulence 2002 !-- is not correlated in space and time, just white noise. This saves 2003 !-- computations power. 2004 ELSE 2005 nux(k) = 0.0_wp 2006 nuy(k) = 0.0_wp 2007 nuz(k) = 0.0_wp 2008 nvx(k) = 0.0_wp 2009 nvy(k) = 0.0_wp 2010 nvz(k) = 0.0_wp 2011 nwx(k) = 0.0_wp 2012 nwy(k) = 0.0_wp 2013 nwz(k) = 0.0_wp 2014 2015 tu(k) = 0.0_wp 2016 tv(k) = 0.0_wp 2017 tw(k) = 0.0_wp 2018 ENDIF 2175 ! 2176 !-- Determine blending parameter. Within the boundary layer length scales 2177 !-- are constant, while above lengths scales approach gradully zero, 2178 !-- i.e. imposed turbulence is not correlated in space and time, 2179 !-- just white noise, which saves computations power as the loops for the 2180 !-- computation of the filter functions depend on the length scales. 2181 !-- The value decreases to 1/10 after one length scale above the 2182 !-- ABL top. 2183 blend = MIN( 1.0_wp, EXP( d_l * zu(k) - d_l * zi_ribulk ) ) 2184 ! 2185 !-- Assume isotropic turbulence length scales 2186 nux(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2187 nuy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2188 nvx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2189 nvy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2190 nwx(k) = MAX( INT( length_scale * ddx ), 1 ) * blend 2191 nwy(k) = MAX( INT( length_scale * ddy ), 1 ) * blend 2192 ! 2193 !-- Along the vertical direction limit the length scale further by the 2194 !-- boundary-layer depth to assure that no length scales larger than 2195 !-- the boundary-layer depth are used 2196 length_scale_dum = MIN( length_scale, zi_ribulk ) 2197 2198 nuz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend 2199 nvz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend 2200 nwz(k) = MAX( INT( length_scale_dum * ddzw(k) ), 1 ) * blend 2201 ! 2202 !-- Limit time scales, else they become very larger for low wind speed, 2203 !-- imposing long-living inflow perturbations which in turn propagate 2204 !-- further into the model domain. Use u_init and v_init to calculate 2205 !-- the time scales, which will be equal to the inflow profiles, both, 2206 !-- in offline nesting mode or in dirichlet/radiation mode. 2207 tu(k) = MIN( dt_stg_adjust, length_scale / & 2208 ( ABS( u_init(k) ) + 0.1_wp ) ) * blend 2209 tv(k) = MIN( dt_stg_adjust, length_scale / & 2210 ( ABS( v_init(k) ) + 0.1_wp ) ) * blend 2211 ! 2212 !-- Time scale of w-component is a mixture from u- and v-component. 2213 tw(k) = SQRT( tu(k)**2 + tv(k)**2 ) * blend 2214 2019 2215 ENDDO 2020 2216 ! … … 2045 2241 !------------------------------------------------------------------------------! 2046 2242 SUBROUTINE calc_scaling_variables 2047 2048 USE arrays_3d, &2049 ONLY: drho_air2050 2051 USE indices, &2052 ONLY: nx, ny2053 2243 2054 2244 USE surface_mod, &
Note: See TracChangeset
for help on using the changeset viewer.
