Changeset 4102 for palm/trunk/SOURCE/turbulence_closure_mod.f90
 Timestamp:
 Jul 17, 2019 4:00:03 PM (2 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

palm/trunk/SOURCE/turbulence_closure_mod.f90
r4048 r4102 25 25 !  26 26 ! $Id$ 27 !  Modularize setting of boundary conditions for TKE and dissipation 28 !  Neumann boundary condition for TKE at model top is set also in child domain 29 !  Revise setting of Neumann boundary conditions at noncyclic lateral 30 ! boundaries 31 !  Bugfix, set Neumann boundary condition for TKE at vertical wall instead of 32 ! an implicit Dirichlet boundary condition which implied a sink of TKE 33 ! at vertical walls 34 ! 35 ! 4048 20190621 21:00:21Z knoop 27 36 ! write out preprocessor directives; remove tailing whitespaces 28 37 ! … … 192 201 !> @todo Check for random disturbances 193 202 !> @note <Enter notes on the module> 194 ! !203 !! 195 204 MODULE turbulence_closure_mod 196 205 197 206 198 USE arrays_3d, 199 ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, 200 e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, 207 USE arrays_3d, & 208 ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, & 209 e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, & 201 210 te_m, tend, u, v, vpt, w 202 211 203 USE basic_constants_and_equations_mod, 212 USE basic_constants_and_equations_mod, & 204 213 ONLY: g, kappa, lv_d_cp, lv_d_rd, rd_d_rv 205 214 206 USE control_parameters, & 207 ONLY: constant_diffusion, dt_3d, e_init, humidity, & 208 initializing_actions, intermediate_timestep_count, & 209 intermediate_timestep_count_max, km_constant, & 210 les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, & 211 pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 212 timestep_scheme, turbulence_closure, & 213 turbulent_inflow, use_upstream_for_tke, vpt_reference, & 215 USE control_parameters, & 216 ONLY: bc_dirichlet_l, & 217 bc_dirichlet_n, & 218 bc_dirichlet_r, & 219 bc_dirichlet_s, & 220 bc_radiation_l, & 221 bc_radiation_n, & 222 bc_radiation_r, & 223 bc_radiation_s, & 224 child_domain, & 225 constant_diffusion, dt_3d, e_init, humidity, & 226 initializing_actions, intermediate_timestep_count, & 227 intermediate_timestep_count_max, km_constant, & 228 les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, & 229 pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 230 timestep_scheme, turbulence_closure, & 231 turbulent_inflow, use_upstream_for_tke, vpt_reference, & 214 232 ws_scheme_sca, current_timestep_number 215 233 216 USE advec_ws, 234 USE advec_ws, & 217 235 ONLY: advec_s_ws 218 236 219 USE advec_s_bc_mod, 237 USE advec_s_bc_mod, & 220 238 ONLY: advec_s_bc 221 239 222 USE advec_s_pw_mod, 240 USE advec_s_pw_mod, & 223 241 ONLY: advec_s_pw 224 242 225 USE advec_s_up_mod, 243 USE advec_s_up_mod, & 226 244 ONLY: advec_s_up 227 245 228 USE cpulog, 246 USE cpulog, & 229 247 ONLY: cpu_log, log_point_s 230 248 231 USE indices, 232 ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, 249 USE indices, & 250 ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 233 251 wall_flags_0 234 252 235 253 USE kinds 236 254 237 USE ocean_mod, 255 USE ocean_mod, & 238 256 ONLY: prho_reference 239 257 240 258 USE pegrid 241 259 242 USE plant_canopy_model_mod, 260 USE plant_canopy_model_mod, & 243 261 ONLY: pcm_tendency 244 262 245 USE statistics, 263 USE statistics, & 246 264 ONLY: hom, hom_sum, statistic_regions 247 265 266 USE surface_mod, & 267 ONLY: bc_h, & 268 bc_v, & 269 get_topography_top_index_ji, & 270 surf_def_h, & 271 surf_def_v, & 272 surf_lsm_h, & 273 surf_lsm_v, & 274 surf_usm_h, & 275 surf_usm_v 248 276 249 277 IMPLICIT NONE … … 271 299 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< nearwall mixing length 272 300 273 301 ! 302 ! Public variables 274 303 PUBLIC c_0, rans_const_c, rans_const_sigma 275 304 … … 277 306 278 307 PRIVATE 279 280 281 PUBLIC & 282 tcm_check_parameters, & 283 tcm_check_data_output, & 284 tcm_define_netcdf_grid, & 285 tcm_init_arrays, & 286 tcm_init, & 287 tcm_actions, & 288 tcm_prognostic_equations, & 289 tcm_swap_timelevel, & 290 tcm_3d_data_averaging, & 291 tcm_data_output_2d, & 292 tcm_data_output_3d, & 308 ! 309 ! Public subroutines 310 PUBLIC & 311 tcm_boundary_conds, & 312 tcm_check_parameters, & 313 tcm_check_data_output, & 314 tcm_define_netcdf_grid, & 315 tcm_init_arrays, & 316 tcm_init, & 317 tcm_actions, & 318 tcm_prognostic_equations, & 319 tcm_swap_timelevel, & 320 tcm_3d_data_averaging, & 321 tcm_data_output_2d, & 322 tcm_data_output_3d, & 293 323 tcm_diffusivities 294 324 295 325 ! 296 326 ! PALM interfaces: 327 ! Boundary conditions for subgrid TKE and dissipation 328 INTERFACE tcm_boundary_conds 329 MODULE PROCEDURE tcm_boundary_conds 330 END INTERFACE tcm_boundary_conds 331 ! 297 332 ! Input parameter checks to be done in check_parameters 298 333 INTERFACE tcm_check_parameters … … 371 406 CONTAINS 372 407 408 !! 409 ! Description: 410 !  411 !> Check parameters routine for turbulence closure module. 412 !! 413 SUBROUTINE tcm_boundary_conds 414 415 USE pmc_interface, & 416 ONLY : rans_mode_parent 417 418 IMPLICIT NONE 419 420 INTEGER(iwp) :: i !< grid index x direction 421 INTEGER(iwp) :: j !< grid index y direction 422 INTEGER(iwp) :: k !< grid index z direction 423 INTEGER(iwp) :: l !< running index boundary type, for up and downwardfacing walls 424 INTEGER(iwp) :: m !< running index surface elements 425 ! 426 ! Boundary conditions for TKE. 427 IF ( .NOT. constant_diffusion ) THEN 428 ! 429 ! In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls. 430 ! Note, only TKE is prognostic in this case and dissipation is only 431 ! a diagnostic quantity. 432 IF ( .NOT. rans_mode ) THEN 433 ! 434 ! Horizontal walls, upward and downwardfacing 435 DO l = 0, 1 436 !$OMP PARALLEL DO PRIVATE( i, j, k ) 437 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 438 !$ACC PRESENT(bc_h, e_p) 439 DO m = 1, bc_h(l)%ns 440 i = bc_h(l)%i(m) 441 j = bc_h(l)%j(m) 442 k = bc_h(l)%k(m) 443 e_p(k+bc_h(l)%koff,j,i) = e_p(k,j,i) 444 ENDDO 445 ENDDO 446 ! 447 ! Vertical walls 448 DO l = 0, 3 449 !$OMP PARALLEL DO PRIVATE( i, j, k ) 450 !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 451 !$ACC PRESENT(bc_v, e_p) 452 DO m = 1, bc_v(l)%ns 453 i = bc_v(l)%i(m) 454 j = bc_v(l)%j(m) 455 k = bc_v(l)%k(m) 456 e_p(k,j+bc_v(l)%joff,i+bc_v(l)%ioff) = e_p(k,j,i) 457 ENDDO 458 ENDDO 459 ! 460 ! In RANS mode, wall function is used as boundary condition for TKE 461 ELSE 462 ! 463 ! Use wall function within constantflux layer 464 ! Note, grid points listed in bc_h are not included in any calculations in RANS mode and 465 ! are therefore not set here. 466 ! 467 ! Upwardfacing surfaces 468 ! Default surfaces 469 DO m = 1, surf_def_h(0)%ns 470 i = surf_def_h(0)%i(m) 471 j = surf_def_h(0)%j(m) 472 k = surf_def_h(0)%k(m) 473 e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2 474 ENDDO 475 ! 476 ! Natural surfaces 477 DO m = 1, surf_lsm_h%ns 478 i = surf_lsm_h%i(m) 479 j = surf_lsm_h%j(m) 480 k = surf_lsm_h%k(m) 481 e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2 482 ENDDO 483 ! 484 ! Urban surfaces 485 DO m = 1, surf_usm_h%ns 486 i = surf_usm_h%i(m) 487 j = surf_usm_h%j(m) 488 k = surf_usm_h%k(m) 489 e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2 490 ENDDO 491 ! 492 ! Vertical surfaces 493 DO l = 0, 3 494 ! 495 ! Default surfaces 496 DO m = 1, surf_def_v(l)%ns 497 i = surf_def_v(l)%i(m) 498 j = surf_def_v(l)%j(m) 499 k = surf_def_v(l)%k(m) 500 e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2 501 ENDDO 502 ! 503 ! Natural surfaces 504 DO m = 1, surf_lsm_v(l)%ns 505 i = surf_lsm_v(l)%i(m) 506 j = surf_lsm_v(l)%j(m) 507 k = surf_lsm_v(l)%k(m) 508 e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2 509 ENDDO 510 ! 511 ! Urban surfaces 512 DO m = 1, surf_usm_v(l)%ns 513 i = surf_usm_v(l)%i(m) 514 j = surf_usm_v(l)%j(m) 515 k = surf_usm_v(l)%k(m) 516 e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2 517 ENDDO 518 ENDDO 519 ENDIF 520 ! 521 ! Set Neumann boundary condition for TKE at model top. Do this also 522 ! in case of a nested run. 523 !$ACC KERNELS PRESENT(e_p) 524 e_p(nzt+1,:,:) = e_p(nzt,:,:) 525 !$ACC END KERNELS 526 ! 527 ! Nesting case: if parent operates in RANS mode and child in LES mode, 528 ! no TKE is transfered. This case, set Neumann conditions at lateral and 529 ! top child boundaries. 530 ! If not ( both either in RANS or in LES mode ), TKE boundary condition 531 ! is treated in the nesting. 532 If ( child_domain ) THEN 533 IF ( rans_mode_parent .AND. .NOT. rans_mode ) THEN 534 535 e_p(nzt+1,:,:) = e_p(nzt,:,:) 536 IF ( bc_dirichlet_l ) e_p(:,:,nxl1) = e_p(:,:,nxl) 537 IF ( bc_dirichlet_r ) e_p(:,:,nxr+1) = e_p(:,:,nxr) 538 IF ( bc_dirichlet_s ) e_p(:,nys1,:) = e_p(:,nys,:) 539 IF ( bc_dirichlet_n ) e_p(:,nyn+1,:) = e_p(:,nyn,:) 540 541 ENDIF 542 ENDIF 543 ! 544 ! At in and outflow boundaries also set Neumann boundary conditions 545 ! for the SGSTKE. An exception is made for the child domain if 546 ! both parent and child operate in RANS mode. This case no 547 ! lateral Neumann boundary conditions will be set but Dirichlet 548 ! conditions will be set in the nesting. 549 IF ( .NOT. child_domain .AND. .NOT. rans_mode_parent .AND. & 550 .NOT. rans_mode ) THEN 551 IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 552 e_p(:,nys1,:) = e_p(:,nys,:) 553 IF ( rans_tke_e ) diss_p(:,nys1,:) = diss_p(:,nys,:) 554 ENDIF 555 IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 556 e_p(:,nyn+1,:) = e_p(:,nyn,:) 557 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 558 ENDIF 559 IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 560 e_p(:,:,nxl1) = e_p(:,:,nxl) 561 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 562 ENDIF 563 IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 564 e_p(:,:,nxr+1) = e_p(:,:,nxr) 565 IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 566 ENDIF 567 ENDIF 568 ENDIF 569 570 ! 571 ! Boundary conditions for TKE dissipation rate in RANS mode. 572 IF ( rans_tke_e ) THEN 573 ! 574 ! Use wall function within constantflux layer 575 ! Upwardfacing surfaces 576 ! Default surfaces 577 DO m = 1, surf_def_h(0)%ns 578 i = surf_def_h(0)%i(m) 579 j = surf_def_h(0)%j(m) 580 k = surf_def_h(0)%k(m) 581 diss_p(k,j,i) = surf_def_h(0)%us(m)**3 & 582 / ( kappa * surf_def_h(0)%z_mo(m) ) 583 ENDDO 584 ! 585 ! Natural surfaces 586 DO m = 1, surf_lsm_h%ns 587 i = surf_lsm_h%i(m) 588 j = surf_lsm_h%j(m) 589 k = surf_lsm_h%k(m) 590 diss_p(k,j,i) = surf_lsm_h%us(m)**3 & 591 / ( kappa * surf_lsm_h%z_mo(m) ) 592 ENDDO 593 ! 594 ! Urban surfaces 595 DO m = 1, surf_usm_h%ns 596 i = surf_usm_h%i(m) 597 j = surf_usm_h%j(m) 598 k = surf_usm_h%k(m) 599 diss_p(k,j,i) = surf_usm_h%us(m)**3 & 600 / ( kappa * surf_usm_h%z_mo(m) ) 601 ENDDO 602 ! 603 ! Vertical surfaces 604 DO l = 0, 3 605 ! 606 ! Default surfaces 607 DO m = 1, surf_def_v(l)%ns 608 i = surf_def_v(l)%i(m) 609 j = surf_def_v(l)%j(m) 610 k = surf_def_v(l)%k(m) 611 diss_p(k,j,i) = surf_def_v(l)%us(m)**3 & 612 / ( kappa * surf_def_v(l)%z_mo(m) ) 613 ENDDO 614 ! 615 ! Natural surfaces 616 DO m = 1, surf_lsm_v(l)%ns 617 i = surf_lsm_v(l)%i(m) 618 j = surf_lsm_v(l)%j(m) 619 k = surf_lsm_v(l)%k(m) 620 diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 & 621 / ( kappa * surf_lsm_v(l)%z_mo(m) ) 622 ENDDO 623 ! 624 ! Urban surfaces 625 DO m = 1, surf_usm_v(l)%ns 626 i = surf_usm_v(l)%i(m) 627 j = surf_usm_v(l)%j(m) 628 k = surf_usm_v(l)%k(m) 629 diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 & 630 / ( kappa * surf_usm_v(l)%z_mo(m) ) 631 ENDDO 632 ENDDO 633 ! 634 ! Limit change of diss to be between 90% and +100%. Also, set an absolute 635 ! minimum value 636 DO i = nxl, nxr 637 DO j = nys, nyn 638 DO k = nzb, nzt+1 639 diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), & 640 2.0_wp * diss(k,j,i) ), & 641 0.1_wp * diss(k,j,i), & 642 0.0001_wp ) 643 ENDDO 644 ENDDO 645 ENDDO 646 647 648 diss_p(nzt+1,:,:) = diss_p(nzt,:,:) 649 650 ENDIF 651 END SUBROUTINE tcm_boundary_conds 652 373 653 !! 374 654 ! Description: … … 926 1206 USE model_1d_mod, & 927 1207 ONLY: e1d, kh1d, km1d 928 929 USE surface_mod, &930 ONLY: get_topography_top_index_ji931 1208 932 1209 IMPLICIT NONE … … 1713 1990 ONLY: phi_m 1714 1991 1715 USE surface_mod, &1716 ONLY : surf_def_h, surf_lsm_h, surf_usm_h1717 1718 1992 IMPLICIT NONE 1719 1993 … … 2418 2692 USE bulk_cloud_model_mod, & 2419 2693 ONLY: bulk_cloud_model 2420 2421 USE surface_mod, &2422 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &2423 surf_usm_v2424 2694 2425 2695 IMPLICIT NONE … … 3122 3392 ONLY: bulk_cloud_model 3123 3393 3124 USE surface_mod, &3125 ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &3126 surf_usm_v3127 3128 3394 IMPLICIT NONE 3129 3395 … … 3767 4033 ONLY: use_sgs_for_particles, wang_kernel 3768 4034 3769 USE surface_mod, &3770 ONLY : bc_h3771 3772 4035 IMPLICIT NONE 3773 4036 … … 3929 4192 3930 4193 ! 3931 ! Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:) 4194 ! Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:). 4195 ! Note, bc cannot be set in tcm_boundary conditions as the dissipation 4196 ! in LES mode is only a diagnostic quantity. 3932 4197 IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. & 3933 4198 wang_kernel .OR. collision_turbulence ) ) THEN … … 3973 4238 USE particle_attributes, & 3974 4239 ONLY: use_sgs_for_particles, wang_kernel 3975 3976 USE surface_mod, &3977 ONLY : bc_h3978 4240 3979 4241 IMPLICIT NONE … … 4062 4324 ! For each surface type determine start and end index (in case of elevated 4063 4325 ! topography several up/downward facing surfaces may exist. 4326 ! Note, bc cannot be set in tcm_boundary conditions as the dissipation 4327 ! in LES mode is only a diagnostic quantity. 4064 4328 surf_s = bc_h(0)%start_index(j,i) 4065 4329 surf_e = bc_h(0)%end_index(j,i) … … 4322 4586 ONLY: phi_m 4323 4587 4324 USE surface_mod, &4325 ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, &4326 surf_usm_h, surf_usm_v4327 4328 4588 INTEGER(iwp) :: i !< loop index 4329 4589 INTEGER(iwp) :: j !< loop index
Note: See TracChangeset
for help on using the changeset viewer.