Changeset 3347 for palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
 Timestamp:
 Oct 15, 2018 2:21:08 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90
r3289 r3347 25 25 !  26 26 ! $Id$ 27 !  Revise structure of init routine 28 !  introduce new parameters to skip STG for some timesteps 29 !  introduce timedependent parametrization of Reynoldsstress tensor 30 !  Bugfix in allocation of mean_inflow_profiles 31 ! 32 ! 3341 20181015 10:31:27Z suehring 27 33 ! Introduce module parameter for number of inflow profiles 28 34 ! … … 122 128 ! Authors: 123 129 !  124 ! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge 130 ! @author Tobias Gronemeier, Matthias Suehring, Atsushi Inagaki, Micha Gryschka, Christoph Knigge 131 ! 125 132 ! 126 133 ! … … 133 140 !> a time scale for each velocity component. The profiles of length and time 134 141 !> scales, mean u, v, w, e and pt, and all components of the Reynolds stress 135 !> tensor are read from file STG_PROFILES. 142 !> tensor can be either read from file STG_PROFILES, or will be parametrized 143 !> within the boundary layer. 136 144 !> 137 145 !> @todo test restart … … 149 157 150 158 USE arrays_3d, & 151 ONLY: mean_inflow_profiles, u, v,w159 ONLY: mean_inflow_profiles, q, pt, u, v, w, zu, zw 152 160 153 161 USE basic_constants_and_equations_mod, & 154 ONLY: pi162 ONLY: g, kappa, pi 155 163 156 164 USE control_parameters, & … … 159 167 160 168 USE cpulog, & 161 ONLY: cpu_log, log_point , log_point_s169 ONLY: cpu_log, log_point 162 170 163 171 USE indices, & … … 170 178 #endif 171 179 180 USE nesting_offl_mod, & 181 ONLY: zi_ribulk 182 172 183 USE pegrid, & 173 184 ONLY: comm1dx, comm1dy, comm2d, ierr, myidx, myidy, pdims … … 184 195 185 196 186 LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main 187 LOGICAL :: use_syn_turb_gen = .FALSE. !< switch to use synthetic turbulence generator 197 LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main 198 LOGICAL :: parametrize_inflow_turbulence = .FALSE. !< flag indicating that inflow turbulence is either read from file (.FALSE.) or if it parametrized 199 LOGICAL :: use_syn_turb_gen = .FALSE. !< switch to use synthetic turbulence generator 188 200 189 201 INTEGER(iwp) :: id_stg_left !< left lateral boundary core id in case of turbulence generator … … 202 214 INTEGER(iwp) :: nzt_y_stg !< upper bound of z coordinate (required for transposing z on PEs along y) 203 215 204 REAL(wp) :: mc_factor !< mass flux correction factor 205 206 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz !< displacement for MPI_GATHERV 207 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz !< receive count for MPI_GATHERV 208 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz !< displacement for MPI_GATHERV 209 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz !< receive count for MPI_GATHERV 210 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux !< length scale of u in x direction (in gp) 211 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy !< length scale of u in y direction (in gp) 212 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz !< length scale of u in z direction (in gp) 213 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx !< length scale of v in x direction (in gp) 214 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy !< length scale of v in y direction (in gp) 215 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz !< length scale of v in z direction (in gp) 216 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx !< length scale of w in x direction (in gp) 217 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy !< length scale of w in y direction (in gp) 218 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz !< length scale of w in z direction (in gp) 219 220 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed !< seed of random number for rngenerator 221 222 REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation 223 REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation 224 REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation 225 REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation 226 REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation 227 REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation 228 REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u 229 REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v 230 REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w 231 232 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux !< filter function for u in x direction 233 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction 234 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction 235 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx !< filter function for v in x direction 236 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction 237 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction 238 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx !< filter function for w in y direction 239 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction 240 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction 241 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz !< velocity seed for u at xz plane 242 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz !< velocity seed for u at xz plane with new random number 243 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz !< velocity seed for u at yz plane 244 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz !< velocity seed for u at yz plane with new random number 245 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz !< velocity seed for v at xz plane 246 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz !< velocity seed for v at xz plane with new random number 247 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz !< velocity seed for v at yz plane 248 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz !< velocity seed for v at yz plane with new random number 249 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz !< velocity seed for w at xz plane 250 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz !< velocity seed for w at xz plane with new random number 251 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz !< velocity seed for w at yz plane 252 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz !< velocity seed for w at yz plane with new random number 253 216 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_xz !< displacement for MPI_GATHERV 217 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_xz !< receive count for MPI_GATHERV 218 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs_yz !< displacement for MPI_GATHERV 219 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count_yz !< receive count for MPI_GATHERV 220 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nux !< length scale of u in x direction (in gp) 221 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy !< length scale of u in y direction (in gp) 222 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz !< length scale of u in z direction (in gp) 223 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvx !< length scale of v in x direction (in gp) 224 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy !< length scale of v in y direction (in gp) 225 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz !< length scale of v in z direction (in gp) 226 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwx !< length scale of w in x direction (in gp) 227 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy !< length scale of w in y direction (in gp) 228 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz !< length scale of w in z direction (in gp) 229 INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed !< seed of random number for rngenerator 230 231 REAL(wp) :: cosmo_ref = 10.0_wp !< height of first vertical grid level in mesoscale model, used for calculation of scaling parameters 232 REAL(wp) :: dt_stg_adjust = 300.0_wp !< time interval for adjusting turbulence statistics 233 REAL(wp) :: dt_stg_call = 5.0_wp !< time interval for calling synthetic turbulence generator 234 REAL(wp) :: mc_factor !< mass flux correction factor 235 REAL(wp) :: scale_l !< scaling parameter used for turbulence parametrization  Obukhov length 236 REAL(wp) :: scale_us !< scaling parameter used for turbulence parametrization  friction velocity 237 REAL(wp) :: scale_wm !< scaling parameter used for turbulence parametrization  momentum scale 238 REAL(wp) :: time_stg_adjust = 0.0_wp !< time counter for adjusting turbulence information 239 REAL(wp) :: time_stg_call = 0.0_wp !< time counter for calling generator 240 241 242 REAL(wp),DIMENSION(:), ALLOCATABLE :: r11 !< Reynolds parameter 243 REAL(wp),DIMENSION(:), ALLOCATABLE :: r21 !< Reynolds parameter 244 REAL(wp),DIMENSION(:), ALLOCATABLE :: r22 !< Reynolds parameter 245 REAL(wp),DIMENSION(:), ALLOCATABLE :: r31 !< Reynolds parameter 246 REAL(wp),DIMENSION(:), ALLOCATABLE :: r32 !< Reynolds parameter 247 REAL(wp),DIMENSION(:), ALLOCATABLE :: r33 !< Reynolds parameter 248 249 REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation 250 REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation 251 REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation 252 REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation 253 REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation 254 REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation 255 REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u 256 REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v 257 REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w 258 259 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bux !< filter function for u in x direction 260 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction 261 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction 262 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvx !< filter function for v in x direction 263 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction 264 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction 265 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwx !< filter function for w in y direction 266 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction 267 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction 268 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_xz !< velocity seed for u at xz plane 269 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_xz !< velocity seed for u at xz plane with new random number 270 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu_yz !< velocity seed for u at yz plane 271 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo_yz !< velocity seed for u at yz plane with new random number 272 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_xz !< velocity seed for v at xz plane 273 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_xz !< velocity seed for v at xz plane with new random number 274 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv_yz !< velocity seed for v at yz plane 275 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo_yz !< velocity seed for v at yz plane with new random number 276 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_xz !< velocity seed for w at xz plane 277 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_xz !< velocity seed for w at xz plane with new random number 278 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw_yz !< velocity seed for w at yz plane 279 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo_yz !< velocity seed for w at yz plane with new random number 280 281 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_xz !< imposed disturbances at north/south boundary 282 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: dist_yz !< imposed disturbances at north/south boundary 254 283 255 284 ! 256 285 ! PALM interfaces: 286 ! Adjust time and lenght scales, Reynolds stress, and filter functions 287 INTERFACE stg_adjust 288 MODULE PROCEDURE stg_adjust 289 END INTERFACE stg_adjust 290 ! 257 291 ! Input parameter checks to be done in check_parameters 258 292 INTERFACE stg_check_parameters … … 319 353 ! 320 354 ! Public interfaces 321 PUBLIC stg_ check_parameters, stg_header, stg_init, stg_main, stg_parin,&322 stg_ wrd_global, stg_rrd_global355 PUBLIC stg_adjust, stg_check_parameters, stg_header, stg_init, stg_main, & 356 stg_parin, stg_rrd_global, stg_wrd_global 323 357 324 358 ! 325 359 ! Public variables 326 PUBLIC id_stg_left, id_stg_north, id_stg_right, id_stg_south, & 327 use_syn_turb_gen 360 PUBLIC dt_stg_call, dt_stg_adjust, id_stg_left, id_stg_north, & 361 id_stg_right, id_stg_south, parametrize_inflow_turbulence, & 362 time_stg_adjust, time_stg_call, use_syn_turb_gen 328 363 329 364 … … 376 411 message_string = 'Using synthetic turbulence generator ' // & 377 412 'requires %initializing_actions = ' // & 378 '"set_constant_profiles" or "read_restart_data"' 413 '"set_constant_profiles" or "read_restart_data"' //& 414 ', if not offline nesting is applied.' 379 415 CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 ) 380 416 ENDIF … … 382 418 IF ( bc_lr /= 'dirichlet/radiation' ) THEN 383 419 message_string = 'Using synthetic turbulence generator ' // & 384 'requires &bc_lr = "dirichlet/radiation"' 420 'requires &bc_lr = "dirichlet/radiation", ' // & 421 'if not offline nesting is applied.' 385 422 CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 ) 386 423 ENDIF 387 424 IF ( bc_ns /= 'cyclic' ) THEN 388 425 message_string = 'Using synthetic turbulence generator ' // & 389 'requires &bc_ns = "cyclic"' 426 'requires &bc_ns = "cyclic", ' // & 427 'if not offline nesting is applied.' 390 428 CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 ) 391 429 ENDIF … … 425 463 WRITE( io, 3 ) 426 464 ENDIF 465 466 IF ( parametrize_inflow_turbulence ) THEN 467 WRITE( io, 4 ) dt_stg_adjust 468 ELSE 469 WRITE( io, 5 ) 470 ENDIF 427 471 428 472 1 FORMAT (//' Synthetic turbulence generator information:'/ & 429 473 ' '/) 430 2 FORMAT (' synthetic turbulence generator switched on') 431 3 FORMAT (' synthetic turbulence generator switched off') 474 2 FORMAT (' synthetic turbulence generator is switched on') 475 3 FORMAT (' synthetic turbulence generator is switched off') 476 4 FORMAT (' imposed turbulence statistics are parametrized and ' \\ & 477 'ajdusted to boundarylayer development each ', F8.2, ' s' ) 478 5 FORMAT (' imposed turbulence is read from file' ) 432 479 433 480 END SUBROUTINE stg_header … … 443 490 444 491 USE arrays_3d, & 445 ONLY: dzw, ddzw, u_init, v_init, zu , zw492 ONLY: dzw, ddzw, u_init, v_init, zu 446 493 447 494 USE control_parameters, & … … 460 507 IMPLICIT NONE 461 508 462 LOGICAL :: file_stg_exist = .FALSE. !< flag indicati onwhether parameter file for Reynolds stress and length scales exist509 LOGICAL :: file_stg_exist = .FALSE. !< flag indicating whether parameter file for Reynolds stress and length scales exist 463 510 464 511 #if defined( __parallel ) … … 484 531 REAL(wp) :: d21 !< vertical interpolation for a21 485 532 REAL(wp) :: d22 !< vertical interpolation for a22 486 REAL(wp) :: dum_exp !< dummy variable used for exponential vertical decrease of turbulent length scales487 533 REAL(wp) :: luy !< length scale for u in y direction 488 534 REAL(wp) :: luz !< length scale for u in z direction … … 494 540 REAL(wp) :: zz !< height 495 541 496 REAL(wp) :: length_scale_surface !< typical length scale497 REAL(wp) :: r_ii_0 !< correction factor498 REAL(wp) :: time_scale !< typical time scale499 REAL(wp) :: length_scale_z !< typical length scale500 501 REAL(wp),DIMENSION(nzb:nzt+1) :: r11 !< Reynolds parameter502 REAL(wp),DIMENSION(nzb:nzt+1) :: r21 !< Reynolds parameter503 REAL(wp),DIMENSION(nzb:nzt+1) :: r22 !< Reynolds parameter504 REAL(wp),DIMENSION(nzb:nzt+1) :: r31 !< Reynolds parameter505 REAL(wp),DIMENSION(nzb:nzt+1) :: r32 !< Reynolds parameter506 REAL(wp),DIMENSION(nzb:nzt+1) :: r33 !< Reynolds parameter507 508 542 509 543 #if defined( __parallel ) … … 512 546 513 547 CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' ) 514 515 IF ( .NOT. ALLOCATED( mean_inflow_profiles ) ) &516 ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) )517 518 ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), &519 a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), &520 nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1), &521 nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1), &522 nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), &523 tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) )524 548 525 549 #if defined( __parallel ) … … 615 639 ENDDO 616 640 CALL RANDOM_SEED(put=seed) 617 641 ! 642 ! Allocate required arrays 643 ! mean_inflow profiles must not be allocated in offline nesting 644 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 645 IF ( .NOT. ALLOCATED( mean_inflow_profiles ) ) & 646 ALLOCATE( mean_inflow_profiles(nzb:nzt+1,1:num_mean_inflow_profiles) ) 647 ENDIF 648 649 ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), & 650 a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), & 651 nux(nzb:nzt+1), nuy(nzb:nzt+1), nuz(nzb:nzt+1), & 652 nvx(nzb:nzt+1), nvy(nzb:nzt+1), nvz(nzb:nzt+1), & 653 nwx(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), & 654 r11(nzb:nzt+1), r21(nzb:nzt+1), r22(nzb:nzt+1), & 655 r31(nzb:nzt+1), r32(nzb:nzt+1), r33(nzb:nzt+1), & 656 tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) ) 657 658 ALLOCATE ( dist_xz(nzb:nzt+1,nxlg:nxrg,3) ) 659 ALLOCATE ( dist_yz(nzb:nzt+1,nysg:nyng,3) ) 660 dist_xz = 0.0_wp 661 dist_yz = 0.0_wp 662 ! 618 663 ! Read inflow profile 619 664 ! Height levels of profiles in input profile is as follows: … … 669 714 670 715 CLOSE( 90 ) 671 716 ! 717 ! Calculate coefficient matrix from Reynolds stress tensor 718 ! (Lund rotation) 719 CALL calc_coeff_matrix 720 ! 721 ! No information about turbulence and its length scales are available. 722 ! Instead, parametrize turbulence which is imposed at the boundaries. 723 ! Set flag which indicates that turbulence is parametrized, which is done 724 ! later when energybalance models are already initialized. This is 725 ! because the STG needs information about surface properties such as 726 ! roughness to build 'realistic' turbulence profiles. 672 727 ELSE 673 728 ! 674 ! Setup default turbulent length scales. From the numerical point of 675 ! view the imposed perturbations should not be immediately dissipated 676 ! by the numerics. The numerical dissipation, however, acts on scales 677 ! up to 8 x the grid spacing. For this reason, set the turbulence 678 ! length scale to 8 time the grid spacing. Further, above the boundary 679 ! layer height, set turbulence lenght scales to zero (equivalent to not 680 ! imposing any perturbations) in order to save computational costs. 681 ! Typical length (time) scales of 100 m (s) should be a good compromise 682 ! between all stratrifications. Nearsurface variances are fixed to 683 ! 0.1 m2/s2, vertical fluxes are one order of magnitude smaller. 684 ! Vertical fluxes 685 length_scale_surface = 100.0_wp 686 r_ii_0 = 0.1_wp 687 time_scale = 100.0_wp 688 DO k = nzb+1, nzt+1 689 dum_exp = MERGE( zu(k) / ( 0.3* zu(nzt) ), & 690 0.0_wp, & 691 zu(k) > 0.3 * zu(nzt) & 692 ) 693 length_scale_z = 8.0_wp * MIN( dx, dy, dzw(k) ) 694 695 ! IF ( zu(k) > zi_richardson ) 696 697 nux(k) = MAX( INT( length_scale_z * ddx ), 1 ) 698 nuy(k) = MAX( INT( length_scale_z * ddy ), 1 ) 699 nuz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 ) 700 nvx(k) = MAX( INT( length_scale_z * ddx ), 1 ) 701 nvy(k) = MAX( INT( length_scale_z * ddy ), 1 ) 702 nvz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 ) 703 nwx(k) = MAX( INT( length_scale_z * ddx ), 1 ) 704 nwy(k) = MAX( INT( length_scale_z * ddy ), 1 ) 705 nwz(k) = MAX( INT( length_scale_z * ddzw(k) ), 1 ) 706 707 r11(k) = r_ii_0 * EXP( dum_exp ) 708 r22(k) = r_ii_0 * EXP( dum_exp ) 709 r33(k) = r_ii_0 * EXP( dum_exp ) 710 711 r21(k) = 0.1_wp * r_ii_0 * EXP( dum_exp ) 712 r31(k) = 0.1_wp * r_ii_0 * EXP( dum_exp ) 713 r32(k) = 0.1_wp * r_ii_0 * EXP( dum_exp ) 714 715 tu(k) = time_scale 716 tv(k) = time_scale 717 tw(k) = time_scale 718 719 ENDDO 720 721 nux(nzb) = nux(nzb+1) 722 nuy(nzb) = nuy(nzb+1) 723 nuz(nzb) = nuz(nzb+1) 724 nvx(nzb) = nvx(nzb+1) 725 nvy(nzb) = nvy(nzb+1) 726 nvz(nzb) = nvz(nzb+1) 727 nwx(nzb) = nwx(nzb+1) 728 nwy(nzb) = nwy(nzb+1) 729 nwz(nzb) = nwz(nzb+1) 730 731 r11(nzb) = r11(nzb+1) 732 r22(nzb) = r22(nzb+1) 733 r33(nzb) = r33(nzb+1) 734 735 r21(nzb) = r11(nzb+1) 736 r31(nzb) = r31(nzb+1) 737 r32(nzb) = r32(nzb+1) 738 739 tu(nzb) = time_scale 740 tv(nzb) = time_scale 741 tw(nzb) = time_scale 742 729 ! Set flag indicating that turbulence is parametrized 730 parametrize_inflow_turbulence = .TRUE. 731 ! 732 ! Determine boundarylayer depth, which is used to initialize lenght 733 ! scales 734 CALL calc_scaling_variables 735 ! 736 ! Initialize lenght and time scales, which in turn are used 737 ! to initialize the filter functions. 738 CALL calc_length_and_time_scale 739 743 740 ENDIF 744 741 745 742 ! 746 ! Assign initial profiles 743 ! Assign initial profiles. Note, this is only required if turbulent 744 ! inflow from the left is desired, not in case of any of the 745 ! nesting (offline or self nesting) approaches. 747 746 IF ( .NOT. nesting_offline .AND. .NOT. child_domain ) THEN 748 747 u_init = mean_inflow_profiles(:,1) … … 751 750 e_init = MAXVAL( mean_inflow_profiles(:,5) ) 752 751 ENDIF 753 ! 754 ! Calculate coefficient matrix from Reynolds stress tensor (Lund rotation) 755 DO k = nzb, nzt+1 756 IF ( r11(k) > 0.0_wp ) THEN 757 a11(k) = SQRT( r11(k) ) 758 a21(k) = r21(k) / a11(k) 759 ELSE 760 a11(k) = 0.0_wp 761 a21(k) = 0.0_wp 762 ENDIF 763 764 a22(k) = r22(k)  a21(k)**2 765 IF ( a22(k) > 0.0_wp ) THEN 766 a22(k) = SQRT( a22(k) ) 767 ELSE 768 a22(k) = 0.0_wp 769 ENDIF 770 771 ! 772 ! a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11, 773 ! d21, d22) because of different vertical grid 774 IF ( k .le. nzt ) THEN 775 d11 = 0.5_wp * ( r11(k) + r11(k+1) ) 776 IF ( d11 > 0.0_wp ) THEN 777 d11 = SQRT( d11 ) 778 d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11 779 a31(k) = r31(k) / d11 780 ELSE 781 d21 = 0.0_wp 782 a31(k) = 0.0_wp 783 ENDIF 784 785 d22 = 0.5_wp * ( r22(k) + r22(k+1) )  d21 ** 2 786 IF ( d22 > 0.0_wp ) THEN 787 a32(k) = ( r32(k)  d21 * a31(k) ) / SQRT( d22 ) 788 ELSE 789 a32(k) = 0.0_wp 790 ENDIF 791 792 a33(k) = r33(k)  a31(k) ** 2  a32(k) ** 2 793 IF ( a33(k) > 0.0_wp ) THEN 794 a33(k) = SQRT( a33(k) ) 795 ELSE 796 a33(k) = 0.0_wp 797 ENDIF 798 ELSE 799 a31(k) = a31(k1) 800 a32(k) = a32(k1) 801 a33(k) = a33(k1) 802 ENDIF 803 804 ENDDO 752 805 753 ! 806 754 ! Define the size of the filter functions and allocate them. … … 961 909 ENDDO 962 910 963 END SUBROUTINE stg_filter_func911 END SUBROUTINE stg_filter_func 964 912 965 913 … … 977 925 978 926 979 NAMELIST /stg_par/ use_syn_turb_gen927 NAMELIST /stg_par/ dt_stg_adjust, dt_stg_call, use_syn_turb_gen 980 928 981 929 line = ' ' … … 1108 1056 REAL(wp) :: volume_flow_l !< local mass flux through lateral boundary 1109 1057 1110 REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg,5) :: dist_xz !< imposed disturbances at north/south boundary1111 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5) :: dist_yz !< imposed disturbances at left/right boundary1112 1113 1114 1058 CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' ) 1115 1059 1116 1060 ! 1117 1061 ! Calculate time step which is needed for filter functions 1118 dt_stg = dt_3d * weight_substep(intermediate_timestep_count) 1119 1062 dt_stg = MAX( dt_3d, dt_stg_call ) !* weight_substep(intermediate_timestep_count) 1120 1063 ! 1121 1064 ! Initial value of fu, fv, fw … … 1145 1088 velocity_seed_initialized = .TRUE. 1146 1089 ENDIF 1090 1147 1091 ! 1148 1092 ! New set of fu, fv, fw … … 1150 1094 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_left ) 1151 1095 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_left ) 1152 1096 1153 1097 IF ( nesting_offline .OR. ( child_domain .AND. rans_mode_parent & 1154 1155 ! 1156 ! 1157 1158 1159 1160 ! 1161 ! 1162 1163 1164 1165 ! 1166 ! 1167 1168 1169 1098 .AND. .NOT. rans_mode ) ) THEN 1099 ! 1100 ! Generate turbulence at right boundary 1101 CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo_yz, id_stg_right ) 1102 CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo_yz, id_stg_right ) 1103 CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo_yz, id_stg_right ) 1104 ! 1105 ! Generate turbulence at north boundary 1106 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_north ) 1107 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_north ) 1108 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_north ) 1109 ! 1110 ! Generate turbulence at south boundary 1111 CALL stg_generate_seed_xz( nux, nuz, bux, buz, fuo_xz, id_stg_south ) 1112 CALL stg_generate_seed_xz( nvx, nvz, bvx, bvz, fvo_xz, id_stg_south ) 1113 CALL stg_generate_seed_xz( nwx, nwz, bwx, bwz, fwo_xz, id_stg_south ) 1170 1114 ENDIF 1115 1171 1116 ! 1172 1117 ! Turbulence generation at left and or right boundary 1173 1118 IF ( myidx == id_stg_left .OR. myidx == id_stg_right ) THEN 1174 1119 1175 1120 DO j = nysg, nyng 1176 1121 DO k = nzb, nzt + 1 … … 1180 1125 fu_yz(k,j) = fuo_yz(k,j) 1181 1126 ELSE 1182 fu_yz(k,j) = fu_yz(k,j) * EXP( pi * dt_stg * 0.5_wp / tu(k) ) + 1127 fu_yz(k,j) = fu_yz(k,j) * EXP( pi * dt_stg * 0.5_wp / tu(k) ) + & 1183 1128 fuo_yz(k,j) * SQRT( 1.0_wp  EXP( pi * dt_stg / tu(k) ) ) 1184 1129 ENDIF … … 1187 1132 fv_yz(k,j) = fvo_yz(k,j) 1188 1133 ELSE 1189 fv_yz(k,j) = fv_yz(k,j) * EXP( pi * dt_stg * 0.5_wp / tv(k) ) + 1134 fv_yz(k,j) = fv_yz(k,j) * EXP( pi * dt_stg * 0.5_wp / tv(k) ) + & 1190 1135 fvo_yz(k,j) * SQRT( 1.0_wp  EXP( pi * dt_stg / tv(k) ) ) 1191 1136 ENDIF … … 1194 1139 fw_yz(k,j) = fwo_yz(k,j) 1195 1140 ELSE 1196 fw_yz(k,j) = fw_yz(k,j) * EXP( pi * dt_stg * 0.5_wp / tw(k) ) + 1141 fw_yz(k,j) = fw_yz(k,j) * EXP( pi * dt_stg * 0.5_wp / tw(k) ) + & 1197 1142 fwo_yz(k,j) * SQRT( 1.0_wp  EXP( pi * dt_stg / tw(k) ) ) 1198 1143 ENDIF … … 1204 1149 dist_yz(k,j,2) = 0.0_wp 1205 1150 dist_yz(k,j,3) = 0.0_wp 1206 ! dist_yz(k,j,4) = 0.0_wp1207 ! dist_yz(k,j,5) = 0.0_wp1208 1151 ELSE 1209 dist_yz(k,j,1) = a11(k) * fu_yz(k,j)1152 dist_yz(k,j,1) = MIN( a11(k) * fu_yz(k,j), 3.0_wp ) 1210 1153 !experimental test of 1.2 1211 dist_yz(k,j,2) = ( SQRT( a22(k) / MAXVAL(a22) ) & 1212 * 1.2_wp ) & 1213 * ( a21(k) * fu_yz(k,j) & 1214 + a22(k) * fv_yz(k,j) ) 1215 dist_yz(k,j,3) = ( SQRT(a33(k) / MAXVAL(a33) ) & 1216 * 1.3_wp ) & 1217 * ( a31(k) * fu_yz(k,j) & 1218 + a32(k) * fv_yz(k,j) & 1219 + a33(k) * fw_yz(k,j) ) 1220 ! Calculation for pt and e not yet implemented 1221 ! dist_yz(k,j,4) = 0.0_wp 1222 ! dist_yz(k,j,5) = 0.0_wp 1154 dist_yz(k,j,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) ) & 1155 * 1.2_wp ) & 1156 * ( a21(k) * fu_yz(k,j) & 1157 + a22(k) * fv_yz(k,j) ), 3.0_wp ) 1158 dist_yz(k,j,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) ) & 1159 * 1.3_wp ) & 1160 * ( a31(k) * fu_yz(k,j) & 1161 + a32(k) * fv_yz(k,j) & 1162 + a33(k) * fw_yz(k,j) ), 3.0_wp ) 1223 1163 ENDIF 1224 1164 … … 1241 1181 1242 1182 #if defined( __parallel ) 1243 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &1244 1, MPI_REAL, MPI_SUM,comm1dy, ierr )1183 CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, 1, MPI_REAL, MPI_SUM, & 1184 comm1dy, ierr ) 1245 1185 #else 1246 1186 mc_factor = mc_factor_l … … 1298 1238 ! Add disturbances 1299 1239 IF ( myidx == id_stg_left ) THEN 1300 1301 1240 DO j = nys, nyn 1302 DO k = nzb , nzt1241 DO k = nzb+1, nzt 1303 1242 u(k,j,0) = ( u(k,j,0) + dist_yz(k,j,1) ) & 1304 1243 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1305 1244 BTEST( wall_flags_0(k,j,0), 1 ) ) 1245 u(k,j,1) = u(k,j,0) 1306 1246 v(k,j,1) = ( v(k,j,1) + dist_yz(k,j,2) ) & 1307 1247 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & … … 1314 1254 ENDIF 1315 1255 IF ( myidx == id_stg_right ) THEN 1316 1317 1256 DO j = nys, nyn 1318 DO k = nzb , nzt1257 DO k = nzb+1, nzt 1319 1258 u(k,j,nxr+1) = ( u(k,j,nxr+1) + dist_yz(k,j,1) ) & 1320 1259 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & … … 1335 1274 ! Turbulence generation at north and south boundary 1336 1275 IF ( myidy == id_stg_north .OR. myidy == id_stg_south ) THEN 1337 1276 1338 1277 DO i = nxlg, nxrg 1339 1278 DO k = nzb, nzt + 1 … … 1351 1290 ELSE 1352 1291 fv_xz(k,i) = fv_xz(k,i) * EXP( pi * dt_stg * 0.5_wp / tv(k) ) + & 1353 1292 fvo_xz(k,i) * SQRT( 1.0_wp  EXP( pi * dt_stg / tv(k) ) ) 1354 1293 ENDIF 1355 1294 1356 1295 IF ( tw(k) == 0.0_wp ) THEN 1357 1296 fw_xz(k,i) = fwo_xz(k,i) 1358 1297 ELSE 1359 1298 fw_xz(k,i) = fw_xz(k,i) * EXP( pi * dt_stg * 0.5_wp / tw(k) ) + & 1360 1299 fwo_xz(k,i) * SQRT( 1.0_wp  EXP( pi * dt_stg / tw(k) ) ) 1361 1300 ENDIF 1362 1301 ! … … 1369 1308 1370 1309 ELSE 1371 dist_xz(k,i,1) = a11(k) * fu_xz(k,i)1310 dist_xz(k,i,1) = MIN( a11(k) * fu_xz(k,i), 3.0_wp ) 1372 1311 !experimental test of 1.2 1373 dist_xz(k,i,2) = ( SQRT( a22(k) / MAXVAL(a22) )&1374 * 1.2_wp )&1375 * ( a21(k) * fu_xz(k,i)&1376 + a22(k) * fv_xz(k,i))1377 dist_xz(k,i,3) = ( SQRT(a33(k) / MAXVAL(a33) )&1378 * 1.3_wp )&1379 * ( a31(k) * fu_xz(k,i)&1380 + a32(k) * fv_xz(k,i)&1381 + a33(k) * fw_xz(k,i))1312 dist_xz(k,i,2) = MIN( ( SQRT( a22(k) / MAXVAL(a22) ) & 1313 * 1.2_wp ) & 1314 * ( a21(k) * fu_xz(k,i) & 1315 + a22(k) * fv_xz(k,i) ), 3.0_wp ) 1316 dist_xz(k,i,3) = MIN( ( SQRT(a33(k) / MAXVAL(a33) ) & 1317 * 1.3_wp ) & 1318 * ( a31(k) * fu_xz(k,i) & 1319 + a32(k) * fv_xz(k,i) & 1320 + a33(k) * fw_xz(k,i) ), 3.0_wp ) 1382 1321 ENDIF 1383 1322 1384 1323 ENDDO 1385 1324 ENDDO 1325 1386 1326 ! 1387 1327 ! Mass flux correction following Kim et al. (2013) … … 1425 1365 1426 1366 DO i = nxl, nxr 1427 DO k = nzb , nzt1367 DO k = nzb+1, nzt 1428 1368 u(k,1,i) = ( u(k,1,i) + dist_xz(k,i,1) ) & 1429 1369 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & … … 1432 1372 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1433 1373 BTEST( wall_flags_0(k,0,i), 2 ) ) 1374 v(k,1,i) = v(k,0,i) 1434 1375 w(k,1,i) = ( w(k,1,i) + dist_xz(k,i,3) ) & 1435 1376 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & … … 1441 1382 1442 1383 DO i = nxl, nxr 1443 DO k = nzb , nzt1444 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) 1445 * mc_factor * MERGE( 1.0_wp, 0.0_wp, 1384 DO k = nzb+1, nzt 1385 u(k,nyn+1,i) = ( u(k,nyn+1,i) + dist_xz(k,i,1) ) & 1386 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1446 1387 BTEST( wall_flags_0(k,nyn+1,i), 1 ) ) 1447 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i, 1) )&1448 * mc_factor * MERGE( 1.0_wp, 0.0_wp, 1388 v(k,nyn+1,i) = ( v(k,nyn+1,i) + dist_xz(k,i,2) ) & 1389 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1449 1390 BTEST( wall_flags_0(k,nyn+1,i), 2 ) ) 1450 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i, 1) )&1451 * mc_factor * MERGE( 1.0_wp, 0.0_wp, 1391 w(k,nyn+1,i) = ( w(k,nyn+1,i) + dist_xz(k,i,3) ) & 1392 * mc_factor * MERGE( 1.0_wp, 0.0_wp, & 1452 1393 BTEST( wall_flags_0(k,nyn+1,i), 3 ) ) 1453 1394 ENDDO … … 1455 1396 ENDIF 1456 1397 ENDIF 1398 ! 1399 ! Finally, set time counter for calling STG to zero 1400 time_stg_call = 0.0_wp 1457 1401 1458 1402 CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' ) … … 1493 1437 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1494 1438 1495 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_y !< filter func in ydir1496 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_z !< filter func in zdir1439 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_y !< filter function in ydirection 1440 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_z !< filter function in zdirection 1497 1441 REAL(wp), DIMENSION(nzb_x_stg:nzt_x_stg+1,nysg:nyng) :: f_n_l !< local velocity seed 1498 1442 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng) :: f_n !< velocity seed … … 1613 1557 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x !< length scale in xdirection 1614 1558 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in zdirection 1615 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2 !< n_ y*21559 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_x2 !< n_x*2 1616 1560 INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2 !< n_z*2 1617 1561 … … 1620 1564 REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number 1621 1565 1622 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_x !< filter func in ydir1623 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_z !< filter func in zdir1566 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_x !< filter function in xdirection 1567 REAL(wp), DIMENSION(merg:merg,nzb:nzt+1) :: b_z !< filter function in zdirection 1624 1568 REAL(wp), DIMENSION(nzb_y_stg:nzt_y_stg+1,nxlg:nxrg) :: f_n_l !< local velocity seed 1625 1569 REAL(wp), DIMENSION(nzb:nzt+1,nxlg:nxrg) :: f_n !< velocity seed … … 1712 1656 END SUBROUTINE stg_generate_seed_xz 1713 1657 1658 !! 1659 ! Description: 1660 !  1661 !> Parametrization of the Reynolds stress tensor, following the parametrization 1662 !> described in Rotach et al. (1996), which is applied in stateoftheart 1663 !> dispserion modelling. Please note, the parametrization does not distinguish 1664 !> between alongwind and crosswind turbulence. 1665 !! 1666 SUBROUTINE parametrize_reynolds_stress 1667 1668 USE arrays_3d, & 1669 ONLY: zu 1670 1671 IMPLICIT NONE 1672 1673 INTEGER(iwp) :: k !< loop index in zdirection 1674 1675 REAL(wp) :: zzi !< ratio of z/zi 1676 1677 ! 1678 ! 1679 DO k = nzb+1, nzt+1 1680 1681 IF ( zu(k) <= zi_ribulk ) THEN 1682 ! 1683 ! Determine normalized height coordinate 1684 zzi = zu(k) / zi_ribulk 1685 ! 1686 ! u'u' and v'v'. Assume isotropy. Note, add a small negative number 1687 ! to the denominator, else the mergefunction can crash if scale_l is 1688 ! zero. 1689 r11(k) = scale_us**2 * ( & 1690 MERGE( 0.35_wp * ( & 1691  zi_ribulk / ( kappa * scale_l  10E4_wp ) & 1692 )**( 2.0_wp / 3.0_wp ), & 1693 0.0_wp, & 1694 scale_l < 0.0_wp ) & 1695 + 5.0_wp  4.0_wp * zzi & 1696 ) 1697 1698 r22(k) = r11(k) 1699 ! 1700 ! w'w' 1701 r33(k) = scale_wm**2 * ( & 1702 1.5_wp * zzi**( 2.0_wp / 3.0_wp ) * EXP( 2.0_wp * zzi ) & 1703 + ( 1.7_wp  zzi ) * ( scale_us / scale_wm )**2 & 1704 ) 1705 ! 1706 ! u'w' and v'w'. Assume isotropy. 1707 r31(k) =  scale_us**2 * ( & 1708 1.0_wp  EXP( 3.0_wp * ( zzi  1.0_wp ) ) & 1709 ) 1710 1711 r32(k) = r31(k) 1712 ! 1713 ! For u'v' no parametrization exist so far  ?. For simplicity assume 1714 ! a similar profile as for u'w'. 1715 r21(k) = r31(k) 1716 ! 1717 ! Above the boundary layer, assmume laminar flow conditions. 1718 ELSE 1719 r11(k) = 10E8_wp 1720 r22(k) = 10E8_wp 1721 r33(k) = 10E8_wp 1722 r21(k) = 10E8_wp 1723 r31(k) = 10E8_wp 1724 r32(k) = 10E8_wp 1725 ENDIF 1726 ! write(9,*) zu(k), r11(k), r33(k), r31(k), zi_ribulk, scale_us, scale_wm, scale_l 1727 ENDDO 1728 1729 ! 1730 ! Set bottom boundary condition 1731 r11(nzb) = r11(nzb+1) 1732 r22(nzb) = r22(nzb+1) 1733 r33(nzb) = r33(nzb+1) 1734 1735 r21(nzb) = r11(nzb+1) 1736 r31(nzb) = r31(nzb+1) 1737 r32(nzb) = r32(nzb+1) 1738 1739 1740 END SUBROUTINE parametrize_reynolds_stress 1741 1742 !! 1743 ! Description: 1744 !  1745 !> Calculate the coefficient matrix from the Lund rotation. 1746 !! 1747 SUBROUTINE calc_coeff_matrix 1748 1749 IMPLICIT NONE 1750 1751 INTEGER(iwp) :: k !< loop index in zdirection 1752 1753 ! 1754 ! Calculate coefficient matrix. Split loops to allow for loop vectorization. 1755 DO k = nzb+1, nzt+1 1756 IF ( r11(k) > 0.0_wp ) THEN 1757 a11(k) = SQRT( r11(k) ) 1758 a21(k) = r21(k) / a11(k) 1759 a31(k) = r31(k) / a11(k) 1760 ELSE 1761 a11(k) = 10E8_wp 1762 a21(k) = 10E8_wp 1763 a31(k) = 10E8_wp 1764 ENDIF 1765 ENDDO 1766 DO k = nzb+1, nzt+1 1767 a22(k) = r22(k)  a21(k)**2 1768 IF ( a22(k) > 0.0_wp ) THEN 1769 a22(k) = SQRT( a22(k) ) 1770 a32(k) = r32(k)  a21(k) * a31(k) / a22(k) 1771 ELSE 1772 a22(k) = 10E8_wp 1773 a32(k) = 10E8_wp 1774 ENDIF 1775 ENDDO 1776 DO k = nzb+1, nzt+1 1777 a33(k) = r33(k)  a31(k)**2  a32(k)**2 1778 IF ( a33(k) > 0.0_wp ) THEN 1779 a33(k) = SQRT( a33(k) ) 1780 ELSE 1781 a33(k) = 10E8_wp 1782 ENDIF 1783 ENDDO 1784 ! 1785 ! Set bottom boundary condition 1786 a11(nzb) = a11(nzb+1) 1787 a22(nzb) = a22(nzb+1) 1788 a21(nzb) = a21(nzb+1) 1789 a33(nzb) = a33(nzb+1) 1790 a31(nzb) = a31(nzb+1) 1791 a32(nzb) = a32(nzb+1) 1792 1793 END SUBROUTINE calc_coeff_matrix 1794 1795 !! 1796 ! Description: 1797 !  1798 !> This routine controls the readjustment of the turbulence statistics used 1799 !> for generating turbulence at the lateral boundaries. 1800 !! 1801 SUBROUTINE stg_adjust 1802 1803 IMPLICIT NONE 1804 1805 INTEGER(iwp) :: k 1806 1807 CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' ) 1808 ! 1809 ! Compute mean boundary layer height according to RichardsonBulk 1810 ! criterion using the inflow profiles. Further velocity scale as well as 1811 ! mean friction velocity are calculated. 1812 CALL calc_scaling_variables 1813 ! 1814 ! Set length and time scales depending on boundarylayer height 1815 CALL calc_length_and_time_scale 1816 ! 1817 ! Parametrize Reynoldsstress tensor, diagonal elements as well 1818 ! as r21 (v'u'), r31 (w'u'), r32 (w'v'). Parametrization follows 1819 ! Rotach et al. (1996) and is based on boundarylayer depth, 1820 ! friction velocity and velocity scale. 1821 CALL parametrize_reynolds_stress 1822 ! 1823 ! Calculate coefficient matrix from Reynolds stress tensor 1824 ! (Lund rotation) 1825 CALL calc_coeff_matrix 1826 ! 1827 ! Determine filter functions on basis of updated length scales 1828 CALL stg_filter_func( nux, bux ) !filter ux 1829 CALL stg_filter_func( nuy, buy ) !filter uy 1830 CALL stg_filter_func( nuz, buz ) !filter uz 1831 CALL stg_filter_func( nvx, bvx ) !filter vx 1832 CALL stg_filter_func( nvy, bvy ) !filter vy 1833 CALL stg_filter_func( nvz, bvz ) !filter vz 1834 CALL stg_filter_func( nwx, bwx ) !filter wx 1835 CALL stg_filter_func( nwy, bwy ) !filter wy 1836 CALL stg_filter_func( nwz, bwz ) !filter wz 1837 ! 1838 ! Reset time counter for controlling next adjustment to zero 1839 time_stg_adjust = 0.0_wp 1840 1841 CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' ) 1842 1843 END SUBROUTINE stg_adjust 1844 1845 1846 !! 1847 ! Description: 1848 !  1849 !> Calculates turbuluent length and time scales if these are not available 1850 !> from measurements. 1851 !! 1852 SUBROUTINE calc_length_and_time_scale 1853 1854 USE arrays_3d, & 1855 ONLY: dzw, ddzw, u_init, v_init, zu, zw 1856 1857 USE grid_variables, & 1858 ONLY: ddx, ddy, dx, dy 1859 1860 USE statistics, & 1861 ONLY: hom 1862 1863 IMPLICIT NONE 1864 1865 1866 INTEGER(iwp) :: k !< loop index in zdirection 1867 1868 REAL(wp) :: length_scale !< typical length scale 1869 1870 ! 1871 ! In initial call the boundarylayer depth can be zero. This case, set 1872 ! minimum value for boundarylayer depth, to setup length scales correctly. 1873 zi_ribulk = MAX( zi_ribulk, zw(nzb+2) ) 1874 ! 1875 ! Setup default turbulent length scales. From the numerical point of 1876 ! view the imposed perturbations should not be immediately dissipated 1877 ! by the numerics. The numerical dissipation, however, acts on scales 1878 ! up to 8 x the grid spacing. For this reason, set the turbulence 1879 ! length scale to 8 time the grid spacing. Further, above the boundary 1880 ! layer height, set turbulence lenght scales to zero (equivalent to not 1881 ! imposing any perturbations) in order to save computational costs. 1882 ! Typical time scales are derived by assuming Taylors's hypothesis, 1883 ! using the length scales and the mean profiles of u and vcomponent. 1884 DO k = nzb+1, nzt+1 1885 1886 length_scale = 8.0_wp * MIN( dx, dy, dzw(k) ) 1887 1888 IF ( zu(k) <= zi_ribulk ) THEN 1889 ! 1890 ! Assume isotropic turbulence length scales 1891 nux(k) = MAX( INT( length_scale * ddx ), 1 ) 1892 nuy(k) = MAX( INT( length_scale * ddy ), 1 ) 1893 nuz(k) = MAX( INT( length_scale * ddzw(k) ), 1 ) 1894 nvx(k) = MAX( INT( length_scale * ddx ), 1 ) 1895 nvy(k) = MAX( INT( length_scale * ddy ), 1 ) 1896 nvz(k) = MAX( INT( length_scale * ddzw(k) ), 1 ) 1897 nwx(k) = MAX( INT( length_scale * ddx ), 1 ) 1898 nwy(k) = MAX( INT( length_scale * ddy ), 1 ) 1899 nwz(k) = MAX( INT( length_scale * ddzw(k) ), 1 ) 1900 ! 1901 ! Limit time scales, else they become very larger for low wind speed, 1902 ! imposing longliving inflow perturbations which in turn propagate 1903 ! further into the model domain. Use u_init and v_init to calculate 1904 ! the time scales, which will be equal to the inflow profiles, both, 1905 ! in offline nesting mode or in dirichlet/radiation mode. 1906 tu(k) = MIN( dt_stg_adjust, length_scale / & 1907 ( ABS( u_init(k) ) + 0.1_wp ) ) 1908 tv(k) = MIN( dt_stg_adjust, length_scale / & 1909 ( ABS( v_init(k) ) + 0.1_wp ) ) 1910 ! 1911 ! Time scale of wcomponent is a mixture from u and vcomponent. 1912 tw(k) = SQRT( tu(k)**2 + tv(k)**2 ) 1913 ! 1914 ! Above the boundary layer length scales are zero, i.e. imposed turbulence 1915 ! is not correlated in space and time, just white noise. This saves 1916 ! computations power. 1917 ELSE 1918 nux(k) = 0.0_wp 1919 nuy(k) = 0.0_wp 1920 nuz(k) = 0.0_wp 1921 nvx(k) = 0.0_wp 1922 nvy(k) = 0.0_wp 1923 nvz(k) = 0.0_wp 1924 nwx(k) = 0.0_wp 1925 nwy(k) = 0.0_wp 1926 nwz(k) = 0.0_wp 1927 1928 tu(k) = 0.0_wp 1929 tv(k) = 0.0_wp 1930 tw(k) = 0.0_wp 1931 ENDIF 1932 ENDDO 1933 ! 1934 ! Set bottom boundary condition for the length and time scales 1935 nux(nzb) = nux(nzb+1) 1936 nuy(nzb) = nuy(nzb+1) 1937 nuz(nzb) = nuz(nzb+1) 1938 nvx(nzb) = nvx(nzb+1) 1939 nvy(nzb) = nvy(nzb+1) 1940 nvz(nzb) = nvz(nzb+1) 1941 nwx(nzb) = nwx(nzb+1) 1942 nwy(nzb) = nwy(nzb+1) 1943 nwz(nzb) = nwz(nzb+1) 1944 1945 tu(nzb) = tu(nzb+1) 1946 tv(nzb) = tv(nzb+1) 1947 tw(nzb) = tw(nzb+1) 1948 1949 1950 END SUBROUTINE calc_length_and_time_scale 1951 1952 !! 1953 ! Description: 1954 !  1955 !> Calculate scaling variables which are used for turbulence parametrization 1956 !> according to Rotach et al. (1996). Scaling variables are: friction velocity, 1957 !> boundarylayer depth, momentum velocity scale, and Obukhov length. 1958 !! 1959 SUBROUTINE calc_scaling_variables 1960 1961 1962 USE control_parameters, & 1963 ONLY: bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r, bc_dirichlet_s, & 1964 humidity, pt_surface 1965 1966 USE indices, & 1967 ONLY: nx, ny 1968 1969 USE surface_mod, & 1970 ONLY: surf_def_h, surf_lsm_h, surf_usm_h 1971 1972 IMPLICIT NONE 1973 1974 INTEGER(iwp) :: i !< loop index in xdirection 1975 INTEGER(iwp) :: j !< loop index in ydirection 1976 INTEGER(iwp) :: k !< loop index in zdirection 1977 INTEGER(iwp) :: k_ref !< index in zdirection for reference height 1978 INTEGER(iwp) :: m !< surface element index 1979 1980 REAL(wp) :: friction_vel_l !< mean friction veloctiy on subdomain 1981 REAL(wp) :: shf_mean !< mean surface sensible heat flux 1982 REAL(wp) :: shf_mean_l !< mean surface sensible heat flux on subdomain 1983 REAL(wp) :: u_int !< ucomponent 1984 REAL(wp) :: v_int !< vcomponent 1985 REAL(wp) :: vpt_surface !< nearsurface virtual potential temperature 1986 REAL(wp) :: w_convective !< convective velocity scale 1987 REAL(wp) :: z0_mean !< mean roughness length 1988 REAL(wp) :: z0_mean_l !< mean roughness length on subdomain 1989 1990 ! 1991 ! Mean friction velocity and velocity scale. Therefore, 1992 ! precalculate mean roughness length and surface sensible heat flux 1993 ! in the model domain, which are further used to estimate friction 1994 ! velocity and velocity scale. Note, for z0 linear averaging is applied, 1995 ! even though this is known to unestimate the effective roughness. 1996 ! This need to be revised later. 1997 z0_mean_l = 0.0_wp 1998 shf_mean_l = 0.0_wp 1999 DO m = 1, surf_def_h(0)%ns 2000 z0_mean_l = z0_mean_l + surf_def_h(0)%z0(m) 2001 shf_mean_l = shf_mean_l + surf_def_h(0)%shf(m) 2002 ENDDO 2003 DO m = 1, surf_lsm_h%ns 2004 z0_mean_l = z0_mean_l + surf_lsm_h%z0(m) 2005 shf_mean_l = shf_mean_l + surf_lsm_h%shf(m) 2006 ENDDO 2007 DO m = 1, surf_usm_h%ns 2008 z0_mean_l = z0_mean_l + surf_usm_h%z0(m) 2009 shf_mean_l = shf_mean_l + surf_usm_h%shf(m) 2010 ENDDO 2011 2012 #if defined( __parallel ) 2013 CALL MPI_ALLREDUCE( z0_mean_l, z0_mean, 1, MPI_REAL, MPI_SUM, & 2014 comm2d, ierr ) 2015 CALL MPI_ALLREDUCE( shf_mean_l, shf_mean, 1, MPI_REAL, MPI_SUM, & 2016 comm2d, ierr ) 2017 #else 2018 z0_mean = z0_mean_l 2019 shf_mean = shf_mean_l 2020 #endif 2021 z0_mean = z0_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2022 shf_mean = shf_mean / REAL( ( nx + 1 ) * ( ny + 1 ), KIND = wp ) 2023 2024 ! 2025 ! Note, Inifor does not use logarithmic interpolation of the 2026 ! velocity components near the ground, so that nearsurface 2027 ! wind speed and thus the friction velocity is overestimated. 2028 ! However, friction velocity is used for turbulence 2029 ! parametrization, so that more physically meaningful values are important. 2030 ! Hence, derive friction velocity from wind speed at a reference height. 2031 ! For a first guess use 20 m, which is in the range of the first 2032 ! COSMO vertical level. 2033 k_ref = MINLOC( ABS( zu  cosmo_ref ), DIM = 1 )  1 2034 2035 friction_vel_l = 0.0_wp 2036 IF ( bc_dirichlet_l .OR. bc_dirichlet_r ) THEN 2037 2038 i = MERGE( 1, nxr + 1, bc_dirichlet_l ) 2039 2040 DO j = nys, nyn 2041 ! 2042 ! Note, in u and v component the imposed perturbations 2043 ! from the STG are already included. Check whether this 2044 ! makes any difference compared to using the puremean 2045 ! inflow profiles. 2046 u_int = MERGE( u(k_ref,j,i+1), u(k_ref,j,i), bc_dirichlet_l ) 2047 v_int = v(k_ref,j,i) 2048 ! 2049 ! Calculate friction velocity and sumup. Therefore, assume 2050 ! neutral condtions. 2051 friction_vel_l = friction_vel_l + kappa * & 2052 SQRT( u_int * u_int + v_int * v_int ) / & 2053 LOG( zu(k_ref) / z0_mean ) 2054 2055 ENDDO 2056 2057 ENDIF 2058 2059 IF ( bc_dirichlet_s .OR. bc_dirichlet_n ) THEN 2060 2061 j = MERGE( 1, nyn + 1, bc_dirichlet_s ) 2062 2063 DO i = nxl, nxr 2064 2065 u_int = u(k_ref,j,i) 2066 v_int = MERGE( v(k_ref,j+1,i), v(k_ref,j,i), bc_dirichlet_s ) 2067 2068 friction_vel_l = friction_vel_l + kappa * & 2069 SQRT( u_int * u_int + v_int * v_int ) / & 2070 LOG( zu(k_ref) / z0_mean ) 2071 2072 ENDDO 2073 2074 ENDIF 2075 2076 #if defined( __parallel ) 2077 CALL MPI_ALLREDUCE( friction_vel_l, scale_us, 1, MPI_REAL, MPI_SUM, & 2078 comm2d, ierr ) 2079 #else 2080 scale_us = friction_vel_l 2081 #endif 2082 scale_us = scale_us / REAL( 2 * nx + 2 * ny, KIND = wp ) 2083 2084 ! 2085 ! Compute mean Obukhov length 2086 IF ( shf_mean > 0.0_wp ) THEN 2087 scale_l =  pt_surface / ( kappa * g ) * scale_us**3 / shf_mean 2088 ELSE 2089 scale_l = 0.0_wp 2090 ENDIF 2091 2092 ! 2093 ! Compute mean convective velocity scale. Note, in case the mean heat flux 2094 ! is negative, set convective velocity scale to zero. 2095 IF ( shf_mean > 0.0_wp ) THEN 2096 w_convective = ( g * shf_mean * zi_ribulk / pt_surface )**( 1.0_wp / 3.0_wp ) 2097 ELSE 2098 w_convective = 0.0_wp 2099 ENDIF 2100 ! 2101 ! Finally, in order to consider also neutral or stable stratification, 2102 ! compute momentum velocity scale from u* and convective velocity scale, 2103 ! according to Rotach et al. (1996). 2104 scale_wm = ( scale_us**3 + 0.6_wp * w_convective**3 )**( 1.0_wp / 3.0_wp ) 2105 2106 END SUBROUTINE calc_scaling_variables 2107 1714 2108 END MODULE synthetic_turbulence_generator_mod
Note: See TracChangeset
for help on using the changeset viewer.