[2353]  1  !> @file turbulence_closure_mod.f90 

[2761]  2  !! 

 3  ! This file is part of the PALM model system. 

[2353]  4  ! 

[2761]  5  ! PALM is free software: you can redistribute it and/or modify it under the 

 6  ! terms of the GNU General Public License as published by the Free Software 

 7  ! Foundation, either version 3 of the License, or (at your option) any later 

[2353]  8  ! version. 

 9  ! 

[2761]  10  ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 

[2353]  11  ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

 12  ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

 13  ! 

 14  ! You should have received a copy of the GNU General Public License along with 

 15  ! PALM. If not, see <http://www.gnu.org/licenses/>. 

 16  ! 

[2761]  17  ! Copyright 20172018 Leibniz Universitaet Hannover 

[2353]  18  !! 

 19  ! 

 20  ! Current revisions: 

 21  !  

[2918]  22  ! 

[4110]  23  ! 

[2918]  24  ! Former revisions: 

 25  !  

 26  ! $Id: turbulence_closure_mod.f90 4180 20190821 14:37:54Z scharf $ 

[4177]  27  ! add comment 

 28  ! 

 29  ! 4170 20190819 17:12:31Z gronemeier 

[4170]  30  !  add performance optimizations according to K. Ketelsen 

 31  ! to diffusion_e and tcm_diffusivities_default 

 32  !  bugfix in calculating l_wall for vertical walls 

 33  !  bugfix in using l_wall in initialization (consider wall_adjustment_factor) 

 34  !  always initialize diss and save the dissipation to that array 

 35  ! 

 36  ! 4168 20190816 13:50:17Z suehring 

[4168]  37  ! Replace function get_topography_top_index by topo_top_ind 

 38  ! 

 39  ! 4110 20190722 17:05:21Z suehring 

[4110]  40  ! pass integer flag array as well as boundary flags to WS scalar advection 

 41  ! routine 

 42  ! 

 43  ! 4109 20190722 17:00:34Z suehring 

[4102]  44  !  Modularize setting of boundary conditions for TKE and dissipation 

 45  !  Neumann boundary condition for TKE at model top is set also in child domain 

 46  !  Revise setting of Neumann boundary conditions at noncyclic lateral 

 47  ! boundaries 

 48  !  Bugfix, set Neumann boundary condition for TKE at vertical wall instead of 

 49  ! an implicit Dirichlet boundary condition which implied a sink of TKE 

 50  ! at vertical walls 

 51  ! 

 52  ! 4048 20190621 21:00:21Z knoop 

[3776]  53  ! write out preprocessor directives; remove tailing whitespaces 

 54  ! 

 55  ! 3775 20190304 12:40:20Z gronemeier 

[3775]  56  ! removed unused variables 

[3776]  57  ! 

[3775]  58  ! 3724 20190206 16:28:23Z kanani 

[3724]  59  ! Correct doubleused log_point_s units 

[3776]  60  ! 

[3724]  61  ! 3719 20190206 13:10:18Z kanani 

[3776]  62  ! Changed log_point to log_point_s, otherwise this overlaps with 

[3719]  63  ! 'all progn.equations' cpu measurement. 

[3776]  64  ! 

[3719]  65  ! 3684 20190120 20:20:58Z knoop 

[3646]  66  ! Remove unused variable simulated_time 

[3776]  67  ! 

 68  ! 

[2353]  69  ! Description: 

 70  !  

 71  !> This module contains the available turbulence closures for PALM. 

 72  !> 

 73  !> 

 74  !> @todo test initialization for all possibilities 

[4170]  75  !> @todo add OpenMP directives whereever possible 

[2938]  76  !> @todo Check for random disturbances 

[2353]  77  !> @note <Enter notes on the module> 

[4102]  78  !! 

[2353]  79  MODULE turbulence_closure_mod 

 80  

[3776]  81  

[4102]  82  USE arrays_3d, & 

 83  ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, & 

 84  e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, & 

[2680]  85  te_m, tend, u, v, vpt, w 

[2353]  86  

[4102]  87  USE basic_constants_and_equations_mod, & 

[3361]  88  ONLY: g, kappa, lv_d_cp, lv_d_rd, rd_d_rv 

[3274]  89  

[4102]  90  USE control_parameters, & 

 91  ONLY: bc_dirichlet_l, & 

 92  bc_dirichlet_n, & 

 93  bc_dirichlet_r, & 

 94  bc_dirichlet_s, & 

 95  bc_radiation_l, & 

 96  bc_radiation_n, & 

 97  bc_radiation_r, & 

 98  bc_radiation_s, & 

 99  child_domain, & 

 100  constant_diffusion, dt_3d, e_init, humidity, & 

 101  initializing_actions, intermediate_timestep_count, & 

 102  intermediate_timestep_count_max, km_constant, & 

 103  les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, & 

 104  pt_reference, rans_mode, rans_tke_e, rans_tke_l, & 

 105  timestep_scheme, turbulence_closure, & 

 106  turbulent_inflow, use_upstream_for_tke, vpt_reference, & 

[3430]  107  ws_scheme_sca, current_timestep_number 

[2353]  108  

[4102]  109  USE advec_ws, & 

[2353]  110  ONLY: advec_s_ws 

 111  

[4102]  112  USE advec_s_bc_mod, & 

[2353]  113  ONLY: advec_s_bc 

 114  

[4102]  115  USE advec_s_pw_mod, & 

[2353]  116  ONLY: advec_s_pw 

 117  

[4102]  118  USE advec_s_up_mod, & 

[2353]  119  ONLY: advec_s_up 

 120  

[4102]  121  USE cpulog, & 

[3719]  122  ONLY: cpu_log, log_point_s 

[2353]  123  

[4102]  124  USE indices, & 

[4109]  125  ONLY: advc_flags_s, & 

 126  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, & 

[4168]  127  topo_top_ind, & 

[2680]  128  wall_flags_0 

[2353]  129  

 130  USE kinds 

 131  

[4102]  132  USE ocean_mod, & 

[3294]  133  ONLY: prho_reference 

 134  

[2353]  135  USE pegrid 

 136  

[4102]  137  USE plant_canopy_model_mod, & 

[2761]  138  ONLY: pcm_tendency 

[2353]  139  

[4102]  140  USE statistics, & 

[2353]  141  ONLY: hom, hom_sum, statistic_regions 

[4102]  142  

 143  USE surface_mod, & 

 144  ONLY: bc_h, & 

 145  bc_v, & 

 146  surf_def_h, & 

 147  surf_def_v, & 

 148  surf_lsm_h, & 

 149  surf_lsm_v, & 

 150  surf_usm_h, & 

 151  surf_usm_v 

[2353]  152  

 153  IMPLICIT NONE 

 154  

 155  

[3083]  156  REAL(wp) :: c_0 !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES) 

 157  REAL(wp) :: c_1 !< model constant for RANS mode 

 158  REAL(wp) :: c_2 !< model constant for RANS mode 

[3398]  159  REAL(wp) :: c_3 !< model constant for RANS mode 

[3083]  160  REAL(wp) :: c_4 !< model constant for RANS mode 

 161  REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length 

 162  REAL(wp) :: dsig_e = 1.0_wp !< factor to calculate Ke from Km (1/sigma_e) 

 163  REAL(wp) :: dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss) 

[2353]  164  

[3083]  165  REAL(wp), DIMENSION(0:4) :: rans_const_c = & !< model constants for RANS mode (namelist param) 

[3398]  166  (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standardtkee closure 

[3083]  167  

 168  REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param) 

[3086]  169  (/ 1.0_wp, 1.30_wp /) 

[3083]  170  

[2913]  171  REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length according to Blackadar 

[3182]  172  

[3776]  173  REAL(wp), DIMENSION(:), ALLOCATABLE :: l_grid !< geometric mean of grid sizes dx, dy, dz 

[2353]  174  

[2913]  175  REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< nearwall mixing length 

 176  

[4102]  177  ! 

 178  ! Public variables 

[3083]  179  PUBLIC c_0, rans_const_c, rans_const_sigma 

[2358]  180  

[4048]  181  SAVE 

 182  

 183  PRIVATE 

[4102]  184  ! 

 185  ! Public subroutines 

 186  PUBLIC & 

 187  tcm_boundary_conds, & 

 188  tcm_check_parameters, & 

 189  tcm_check_data_output, & 

 190  tcm_define_netcdf_grid, & 

 191  tcm_init_arrays, & 

 192  tcm_init, & 

 193  tcm_actions, & 

 194  tcm_prognostic_equations, & 

 195  tcm_swap_timelevel, & 

 196  tcm_3d_data_averaging, & 

 197  tcm_data_output_2d, & 

 198  tcm_data_output_3d, & 

[4048]  199  tcm_diffusivities 

 200  

[2353]  201  ! 

[2680]  202  ! PALM interfaces: 

[4102]  203  ! Boundary conditions for subgrid TKE and dissipation 

 204  INTERFACE tcm_boundary_conds 

 205  MODULE PROCEDURE tcm_boundary_conds 

 206  END INTERFACE tcm_boundary_conds 

 207  ! 

[2680]  208  ! Input parameter checks to be done in check_parameters 

 209  INTERFACE tcm_check_parameters 

 210  MODULE PROCEDURE tcm_check_parameters 

 211  END INTERFACE tcm_check_parameters 

[2353]  212  

 213  ! 

 214  ! Data output checks for 2D/3D data to be done in check_parameters 

 215  INTERFACE tcm_check_data_output 

 216  MODULE PROCEDURE tcm_check_data_output 

 217  END INTERFACE tcm_check_data_output 

[2680]  218  

[2353]  219  ! 

[2680]  220  ! Definition of data output quantities 

 221  INTERFACE tcm_define_netcdf_grid 

 222  MODULE PROCEDURE tcm_define_netcdf_grid 

 223  END INTERFACE tcm_define_netcdf_grid 

[2353]  224  

 225  ! 

[4048]  226  ! Initialization of arrays 

 227  INTERFACE tcm_init_arrays 

 228  MODULE PROCEDURE tcm_init_arrays 

 229  END INTERFACE tcm_init_arrays 

[2353]  230  

 231  ! 

[3776]  232  ! Initialization actions 

[2353]  233  INTERFACE tcm_init 

 234  MODULE PROCEDURE tcm_init 

 235  END INTERFACE tcm_init 

[2680]  236  

[2353]  237  ! 

[4048]  238  ! Location specific actions 

 239  INTERFACE tcm_actions 

 240  MODULE PROCEDURE tcm_actions 

 241  MODULE PROCEDURE tcm_actions_ij 

 242  END INTERFACE tcm_actions 

[2353]  243  

 244  ! 

[2680]  245  ! Prognostic equations for TKE and TKE dissipation rate 

[3386]  246  INTERFACE tcm_prognostic_equations 

 247  MODULE PROCEDURE tcm_prognostic_equations 

 248  MODULE PROCEDURE tcm_prognostic_equations_ij 

 249  END INTERFACE tcm_prognostic_equations 

[2353]  250  

[2680]  251  ! 

[4048]  252  ! Swapping of time levels (required for prognostic variables) 

 253  INTERFACE tcm_swap_timelevel 

 254  MODULE PROCEDURE tcm_swap_timelevel 

 255  END INTERFACE tcm_swap_timelevel 

[2353]  256  

[2680]  257  ! 

[4048]  258  ! Averaging of 3D data for output 

 259  INTERFACE tcm_3d_data_averaging 

 260  MODULE PROCEDURE tcm_3d_data_averaging 

 261  END INTERFACE tcm_3d_data_averaging 

[2353]  262  

[2680]  263  ! 

[4048]  264  ! Data output of 2D quantities 

 265  INTERFACE tcm_data_output_2d 

 266  MODULE PROCEDURE tcm_data_output_2d 

 267  END INTERFACE tcm_data_output_2d 

[2353]  268  

[2680]  269  ! 

[4048]  270  ! Data output of 3D data 

 271  INTERFACE tcm_data_output_3d 

 272  MODULE PROCEDURE tcm_data_output_3d 

 273  END INTERFACE tcm_data_output_3d 

[2353]  274  

 275  ! 

[3120]  276  ! Call tcm_diffusivities_default and tcm_diffusivities_dynamic 

[2680]  277  INTERFACE tcm_diffusivities 

 278  MODULE PROCEDURE tcm_diffusivities 

 279  END INTERFACE tcm_diffusivities 

[2353]  280  

[3120]  281  

[2680]  282  CONTAINS 

[2353]  283  

 284  !! 

 285  ! Description: 

 286  !  

[2680]  287  !> Check parameters routine for turbulence closure module. 

[2353]  288  !! 

[4102]  289  SUBROUTINE tcm_boundary_conds 

 290  

 291  USE pmc_interface, & 

 292  ONLY : rans_mode_parent 

 293  

 294  IMPLICIT NONE 

 295  

 296  INTEGER(iwp) :: i !< grid index x direction 

 297  INTEGER(iwp) :: j !< grid index y direction 

 298  INTEGER(iwp) :: k !< grid index z direction 

 299  INTEGER(iwp) :: l !< running index boundary type, for up and downwardfacing walls 

 300  INTEGER(iwp) :: m !< running index surface elements 

 301  ! 

 302  ! Boundary conditions for TKE. 

 303  IF ( .NOT. constant_diffusion ) THEN 

 304  ! 

 305  ! In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls. 

 306  ! Note, only TKE is prognostic in this case and dissipation is only 

 307  ! a diagnostic quantity. 

 308  IF ( .NOT. rans_mode ) THEN 

 309  ! 

 310  ! Horizontal walls, upward and downwardfacing 

 311  DO l = 0, 1 

 312  !$OMP PARALLEL DO PRIVATE( i, j, k ) 

 313  !$ACC PARALLEL LOOP PRIVATE(i, j, k) & 

 314  !$ACC PRESENT(bc_h, e_p) 

 315  DO m = 1, bc_h(l)%ns 

 316  i = bc_h(l)%i(m) 

 317  j = bc_h(l)%j(m) 

 318  k = bc_h(l)%k(m) 

 319  e_p(k+bc_h(l)%koff,j,i) = e_p(k,j,i) 

 320  ENDDO 

 321  ENDDO 

 322  ! 

 323  ! Vertical walls 

 324  DO l = 0, 3 

[4105]  325  ! 

 326  ! Note concerning missing ACC directive for this loop: Even though 

 327  ! the data structure bc_v is present, it may not contain any 

 328  ! allocated arrays in the flat but also in a topography case, 

 329  ! leading to a runtime error. Therefore, omit ACC directives 

 330  ! for this loop, in contrast to the bc_h loop. 

[4102]  331  !$OMP PARALLEL DO PRIVATE( i, j, k ) 

 332  DO m = 1, bc_v(l)%ns 

[4105]  333  i = bc_v(l)%i(m) 

[4102]  334  j = bc_v(l)%j(m) 

 335  k = bc_v(l)%k(m) 

 336  e_p(k,j+bc_v(l)%joff,i+bc_v(l)%ioff) = e_p(k,j,i) 

 337  ENDDO 

 338  ENDDO 

 339  ! 

 340  ! In RANS mode, wall function is used as boundary condition for TKE 

 341  ELSE 

 342  ! 

 343  ! Use wall function within constantflux layer 

 344  ! Note, grid points listed in bc_h are not included in any calculations in RANS mode and 

 345  ! are therefore not set here. 

 346  ! 

 347  ! Upwardfacing surfaces 

 348  ! Default surfaces 

 349  DO m = 1, surf_def_h(0)%ns 

 350  i = surf_def_h(0)%i(m) 

 351  j = surf_def_h(0)%j(m) 

 352  k = surf_def_h(0)%k(m) 

 353  e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2 

 354  ENDDO 

 355  ! 

 356  ! Natural surfaces 

 357  DO m = 1, surf_lsm_h%ns 

 358  i = surf_lsm_h%i(m) 

 359  j = surf_lsm_h%j(m) 

 360  k = surf_lsm_h%k(m) 

 361  e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2 

 362  ENDDO 

 363  ! 

 364  ! Urban surfaces 

 365  DO m = 1, surf_usm_h%ns 

 366  i = surf_usm_h%i(m) 

 367  j = surf_usm_h%j(m) 

 368  k = surf_usm_h%k(m) 

 369  e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2 

 370  ENDDO 

 371  ! 

 372  ! Vertical surfaces 

 373  DO l = 0, 3 

 374  ! 

 375  ! Default surfaces 

 376  DO m = 1, surf_def_v(l)%ns 

 377  i = surf_def_v(l)%i(m) 

 378  j = surf_def_v(l)%j(m) 

 379  k = surf_def_v(l)%k(m) 

 380  e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2 

 381  ENDDO 

 382  ! 

 383  ! Natural surfaces 

 384  DO m = 1, surf_lsm_v(l)%ns 

 385  i = surf_lsm_v(l)%i(m) 

 386  j = surf_lsm_v(l)%j(m) 

 387  k = surf_lsm_v(l)%k(m) 

 388  e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2 

 389  ENDDO 

 390  ! 

 391  ! Urban surfaces 

 392  DO m = 1, surf_usm_v(l)%ns 

 393  i = surf_usm_v(l)%i(m) 

 394  j = surf_usm_v(l)%j(m) 

 395  k = surf_usm_v(l)%k(m) 

 396  e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2 

 397  ENDDO 

 398  ENDDO 

 399  ENDIF 

 400  ! 

 401  ! Set Neumann boundary condition for TKE at model top. Do this also 

 402  ! in case of a nested run. 

 403  !$ACC KERNELS PRESENT(e_p) 

 404  e_p(nzt+1,:,:) = e_p(nzt,:,:) 

 405  !$ACC END KERNELS 

 406  ! 

 407  ! Nesting case: if parent operates in RANS mode and child in LES mode, 

 408  ! no TKE is transfered. This case, set Neumann conditions at lateral and 

 409  ! top child boundaries. 

 410  ! If not ( both either in RANS or in LES mode ), TKE boundary condition 

 411  ! is treated in the nesting. 

 412  If ( child_domain ) THEN 

 413  IF ( rans_mode_parent .AND. .NOT. rans_mode ) THEN 

 414  

 415  e_p(nzt+1,:,:) = e_p(nzt,:,:) 

 416  IF ( bc_dirichlet_l ) e_p(:,:,nxl1) = e_p(:,:,nxl) 

 417  IF ( bc_dirichlet_r ) e_p(:,:,nxr+1) = e_p(:,:,nxr) 

 418  IF ( bc_dirichlet_s ) e_p(:,nys1,:) = e_p(:,nys,:) 

 419  IF ( bc_dirichlet_n ) e_p(:,nyn+1,:) = e_p(:,nyn,:) 

 420  

 421  ENDIF 

 422  ENDIF 

 423  ! 

 424  ! At in and outflow boundaries also set Neumann boundary conditions 

 425  ! for the SGSTKE. An exception is made for the child domain if 

 426  ! both parent and child operate in RANS mode. This case no 

 427  ! lateral Neumann boundary conditions will be set but Dirichlet 

 428  ! conditions will be set in the nesting. 

 429  IF ( .NOT. child_domain .AND. .NOT. rans_mode_parent .AND. & 

 430  .NOT. rans_mode ) THEN 

 431  IF ( bc_dirichlet_s .OR. bc_radiation_s ) THEN 

 432  e_p(:,nys1,:) = e_p(:,nys,:) 

 433  IF ( rans_tke_e ) diss_p(:,nys1,:) = diss_p(:,nys,:) 

 434  ENDIF 

 435  IF ( bc_dirichlet_n .OR. bc_radiation_n ) THEN 

 436  e_p(:,nyn+1,:) = e_p(:,nyn,:) 

 437  IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 

 438  ENDIF 

 439  IF ( bc_dirichlet_l .OR. bc_radiation_l ) THEN 

 440  e_p(:,:,nxl1) = e_p(:,:,nxl) 

 441  IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 

 442  ENDIF 

 443  IF ( bc_dirichlet_r .OR. bc_radiation_r ) THEN 

 444  e_p(:,:,nxr+1) = e_p(:,:,nxr) 

 445  IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 

 446  ENDIF 

 447  ENDIF 

 448  ENDIF 

 449  

 450  ! 

 451  ! Boundary conditions for TKE dissipation rate in RANS mode. 

 452  IF ( rans_tke_e ) THEN 

 453  ! 

 454  ! Use wall function within constantflux layer 

 455  ! Upwardfacing surfaces 

 456  ! Default surfaces 

 457  DO m = 1, surf_def_h(0)%ns 

 458  i = surf_def_h(0)%i(m) 

 459  j = surf_def_h(0)%j(m) 

 460  k = surf_def_h(0)%k(m) 

 461  diss_p(k,j,i) = surf_def_h(0)%us(m)**3 & 

 462  / ( kappa * surf_def_h(0)%z_mo(m) ) 

 463  ENDDO 

 464  ! 

 465  ! Natural surfaces 

 466  DO m = 1, surf_lsm_h%ns 

 467  i = surf_lsm_h%i(m) 

 468  j = surf_lsm_h%j(m) 

 469  k = surf_lsm_h%k(m) 

 470  diss_p(k,j,i) = surf_lsm_h%us(m)**3 & 

 471  / ( kappa * surf_lsm_h%z_mo(m) ) 

 472  ENDDO 

 473  ! 

 474  ! Urban surfaces 

 475  DO m = 1, surf_usm_h%ns 

 476  i = surf_usm_h%i(m) 

 477  j = surf_usm_h%j(m) 

 478  k = surf_usm_h%k(m) 

 479  diss_p(k,j,i) = surf_usm_h%us(m)**3 & 

 480  / ( kappa * surf_usm_h%z_mo(m) ) 

 481  ENDDO 

 482  ! 

 483  ! Vertical surfaces 

 484  DO l = 0, 3 

 485  ! 

 486  ! Default surfaces 

 487  DO m = 1, surf_def_v(l)%ns 

 488  i = surf_def_v(l)%i(m) 

 489  j = surf_def_v(l)%j(m) 

 490  k = surf_def_v(l)%k(m) 

 491  diss_p(k,j,i) = surf_def_v(l)%us(m)**3 & 

 492  / ( kappa * surf_def_v(l)%z_mo(m) ) 

 493  ENDDO 

 494  ! 

 495  ! Natural surfaces 

 496  DO m = 1, surf_lsm_v(l)%ns 

 497  i = surf_lsm_v(l)%i(m) 

 498  j = surf_lsm_v(l)%j(m) 

 499  k = surf_lsm_v(l)%k(m) 

 500  diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 & 

 501  / ( kappa * surf_lsm_v(l)%z_mo(m) ) 

 502  ENDDO 

 503  ! 

 504  ! Urban surfaces 

 505  DO m = 1, surf_usm_v(l)%ns 

 506  i = surf_usm_v(l)%i(m) 

 507  j = surf_usm_v(l)%j(m) 

 508  k = surf_usm_v(l)%k(m) 

 509  diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 & 

 510  / ( kappa * surf_usm_v(l)%z_mo(m) ) 

 511  ENDDO 

 512  ENDDO 

 513  ! 

 514  ! Limit change of diss to be between 90% and +100%. Also, set an absolute 

 515  ! minimum value 

 516  DO i = nxl, nxr 

 517  DO j = nys, nyn 

 518  DO k = nzb, nzt+1 

 519  diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), & 

 520  2.0_wp * diss(k,j,i) ), & 

 521  0.1_wp * diss(k,j,i), & 

 522  0.0001_wp ) 

 523  ENDDO 

 524  ENDDO 

 525  ENDDO 

 526  

 527  diss_p(nzt+1,:,:) = diss_p(nzt,:,:) 

[4170]  528  

[4102]  529  ENDIF 

[4170]  530  

[4102]  531  END SUBROUTINE tcm_boundary_conds 

 532  

 533  !! 

 534  ! Description: 

 535  !  

 536  !> Check parameters routine for turbulence closure module. 

 537  !! 

[2353]  538  SUBROUTINE tcm_check_parameters 

 539  

 540  USE control_parameters, & 

[3241]  541  ONLY: message_string, turbulent_inflow, turbulent_outflow 

[2353]  542  

 543  IMPLICIT NONE 

 544  

 545  ! 

 546  ! Define which turbulence closure is going to be used 

[3545]  547  SELECT CASE ( TRIM( turbulence_closure ) ) 

[2353]  548  

[3545]  549  CASE ( 'dynamic' ) 

 550  les_dynamic = .TRUE. 

 551  

 552  CASE ( 'Moeng_Wyngaard' ) 

 553  les_mw = .TRUE. 

 554  

 555  CASE ( 'TKEl' ) 

 556  rans_tke_l = .TRUE. 

 557  rans_mode = .TRUE. 

 558  

 559  CASE ( 'TKEe' ) 

 560  rans_tke_e = .TRUE. 

 561  rans_mode = .TRUE. 

 562  

 563  CASE DEFAULT 

 564  message_string = 'Unknown turbulence closure: ' // & 

 565  TRIM( turbulence_closure ) 

 566  CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 ) 

 567  

 568  END SELECT 

[3083]  569  ! 

[3545]  570  ! Set variables for RANS mode or LES mode 

 571  IF ( rans_mode ) THEN 

 572  ! 

[3083]  573  ! Assign values to constants for RANS mode 

 574  dsig_e = 1.0_wp / rans_const_sigma(1) 

 575  dsig_diss = 1.0_wp / rans_const_sigma(2) 

[2353]  576  

[3083]  577  c_0 = rans_const_c(0) 

 578  c_1 = rans_const_c(1) 

 579  c_2 = rans_const_c(2) 

[3398]  580  c_3 = rans_const_c(3) !> @todo clarify how to switch between different models 

[3083]  581  c_4 = rans_const_c(4) 

 582  

 583  IF ( turbulent_inflow .OR. turbulent_outflow ) THEN 

 584  message_string = 'turbulent inflow/outflow is not yet '// & 

 585  'implemented for RANS mode' 

 586  CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 ) 

 587  ENDIF 

 588  

[2353]  589  message_string = 'RANS mode is still in development! ' // & 

 590  '&Not all features of PALM are yet compatible '// & 

 591  'with RANS mode. &Use at own risk!' 

[3083]  592  CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 ) 

[2353]  593  

 594  ELSE 

[3545]  595  ! 

 596  ! LES mode 

 597  c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980) 

[2353]  598  

[3776]  599  dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead 

[3083]  600  !of K_e which is used in RANS mode 

 601  

[2353]  602  ENDIF 

 603  

 604  END SUBROUTINE tcm_check_parameters 

 605  

 606  !! 

[2680]  607  ! Description: 

 608  !  

 609  !> Check data output. 

 610  !! 

[3241]  611  SUBROUTINE tcm_check_data_output( var, unit ) 

[3776]  612  

[2680]  613  IMPLICIT NONE 

 614  

[3083]  615  CHARACTER (LEN=*) :: unit !< unit of output variable 

 616  CHARACTER (LEN=*) :: var !< name of output variable 

[2680]  617  

 618  

 619  SELECT CASE ( TRIM( var ) ) 

 620  

 621  CASE ( 'diss' ) 

 622  unit = 'm2/s3' 

 623  

 624  CASE ( 'kh', 'km' ) 

 625  unit = 'm2/s' 

 626  

 627  CASE DEFAULT 

 628  unit = 'illegal' 

 629  

 630  END SELECT 

 631  

 632  END SUBROUTINE tcm_check_data_output 

 633  

 634  

 635  !! 

 636  ! Description: 

 637  !  

 638  !> Define appropriate grid for netcdf variables. 

 639  !> It is called out from subroutine netcdf. 

 640  !! 

 641  SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z ) 

[3776]  642  

[2680]  643  IMPLICIT NONE 

 644  

[3083]  645  CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable 

 646  CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable 

 647  CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable 

 648  CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable 

 649  

 650  LOGICAL, INTENT(OUT) :: found !< flag if output variable is found 

 651  

[2680]  652  found = .TRUE. 

 653  

[2353]  654  ! 

[2680]  655  ! Check for the grid 

 656  SELECT CASE ( TRIM( var ) ) 

 657  

 658  CASE ( 'diss', 'diss_xy', 'diss_xz', 'diss_yz' ) 

 659  grid_x = 'x' 

 660  grid_y = 'y' 

 661  grid_z = 'zu' 

 662  

 663  CASE ( 'kh', 'kh_xy', 'kh_xz', 'kh_yz' ) 

 664  grid_x = 'x' 

 665  grid_y = 'y' 

 666  grid_z = 'zu' 

 667  

 668  CASE ( 'km', 'km_xy', 'km_xz', 'km_yz' ) 

 669  grid_x = 'x' 

 670  grid_y = 'y' 

 671  grid_z = 'zu' 

 672  

 673  CASE DEFAULT 

 674  found = .FALSE. 

 675  grid_x = 'none' 

 676  grid_y = 'none' 

 677  grid_z = 'none' 

 678  

 679  END SELECT 

 680  

 681  END SUBROUTINE tcm_define_netcdf_grid 

 682  

 683  

 684  !! 

[2353]  685  ! Description: 

 686  !  

[2680]  687  !> Average 3D data. 

[2353]  688  !! 

 689  SUBROUTINE tcm_3d_data_averaging( mode, variable ) 

 690  

[3776]  691  

[2353]  692  USE averaging, & 

[2680]  693  ONLY: diss_av, kh_av, km_av 

[2353]  694  

[2680]  695  USE control_parameters, & 

 696  ONLY: average_count_3d 

[2353]  697  

 698  IMPLICIT NONE 

 699  

[3083]  700  CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average' 

 701  CHARACTER (LEN=*) :: variable !< name of variable 

[2353]  702  

[3083]  703  INTEGER(iwp) :: i !< loop index 

 704  INTEGER(iwp) :: j !< loop index 

 705  INTEGER(iwp) :: k !< loop index 

[2353]  706  

 707  IF ( mode == 'allocate' ) THEN 

 708  

 709  SELECT CASE ( TRIM( variable ) ) 

 710  

 711  CASE ( 'diss' ) 

 712  IF ( .NOT. ALLOCATED( diss_av ) ) THEN 

[2680]  713  ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

[2353]  714  ENDIF 

 715  diss_av = 0.0_wp 

 716  

[2680]  717  CASE ( 'kh' ) 

 718  IF ( .NOT. ALLOCATED( kh_av ) ) THEN 

 719  ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 720  ENDIF 

 721  kh_av = 0.0_wp 

 722  

 723  CASE ( 'km' ) 

 724  IF ( .NOT. ALLOCATED( km_av ) ) THEN 

 725  ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 726  ENDIF 

 727  km_av = 0.0_wp 

 728  

[2353]  729  CASE DEFAULT 

 730  CONTINUE 

 731  

 732  END SELECT 

 733  

 734  ELSEIF ( mode == 'sum' ) THEN 

 735  

 736  SELECT CASE ( TRIM( variable ) ) 

 737  

 738  CASE ( 'diss' ) 

[3776]  739  IF ( ALLOCATED( diss_av ) ) THEN 

[3004]  740  DO i = nxlg, nxrg 

 741  DO j = nysg, nyng 

 742  DO k = nzb, nzt+1 

 743  diss_av(k,j,i) = diss_av(k,j,i) + diss(k,j,i) 

 744  ENDDO 

[2353]  745  ENDDO 

 746  ENDDO 

[3004]  747  ENDIF 

[2353]  748  

[2680]  749  CASE ( 'kh' ) 

[3004]  750  IF ( ALLOCATED( kh_av ) ) THEN 

 751  DO i = nxlg, nxrg 

 752  DO j = nysg, nyng 

 753  DO k = nzb, nzt+1 

 754  kh_av(k,j,i) = kh_av(k,j,i) + kh(k,j,i) 

 755  ENDDO 

[2680]  756  ENDDO 

 757  ENDDO 

[3004]  758  ENDIF 

[2680]  759  

 760  CASE ( 'km' ) 

[3004]  761  IF ( ALLOCATED( km_av ) ) THEN 

 762  DO i = nxlg, nxrg 

 763  DO j = nysg, nyng 

 764  DO k = nzb, nzt+1 

 765  km_av(k,j,i) = km_av(k,j,i) + km(k,j,i) 

 766  ENDDO 

[2680]  767  ENDDO 

 768  ENDDO 

[3004]  769  ENDIF 

[2680]  770  

[2353]  771  CASE DEFAULT 

 772  CONTINUE 

 773  

 774  END SELECT 

 775  

 776  ELSEIF ( mode == 'average' ) THEN 

 777  

 778  SELECT CASE ( TRIM( variable ) ) 

 779  

 780  CASE ( 'diss' ) 

[3004]  781  IF ( ALLOCATED( diss_av ) ) THEN 

 782  DO i = nxlg, nxrg 

 783  DO j = nysg, nyng 

 784  DO k = nzb, nzt+1 

[3776]  785  diss_av(k,j,i) = diss_av(k,j,i) & 

[3004]  786  / REAL( average_count_3d, KIND=wp ) 

 787  ENDDO 

[2353]  788  ENDDO 

 789  ENDDO 

[3004]  790  ENDIF 

[2353]  791  

[2680]  792  CASE ( 'kh' ) 

[3004]  793  IF ( ALLOCATED( kh_av ) ) THEN 

 794  DO i = nxlg, nxrg 

 795  DO j = nysg, nyng 

 796  DO k = nzb, nzt+1 

[3776]  797  kh_av(k,j,i) = kh_av(k,j,i) & 

[3004]  798  / REAL( average_count_3d, KIND=wp ) 

 799  ENDDO 

[2680]  800  ENDDO 

 801  ENDDO 

[3004]  802  ENDIF 

[2680]  803  

 804  CASE ( 'km' ) 

[3004]  805  IF ( ALLOCATED( km_av ) ) THEN 

 806  DO i = nxlg, nxrg 

 807  DO j = nysg, nyng 

 808  DO k = nzb, nzt+1 

[3776]  809  km_av(k,j,i) = km_av(k,j,i) & 

[3004]  810  / REAL( average_count_3d, KIND=wp ) 

 811  ENDDO 

[2680]  812  ENDDO 

 813  ENDDO 

[3004]  814  ENDIF 

[2680]  815  

[2353]  816  END SELECT 

 817  

 818  ENDIF 

 819  

 820  END SUBROUTINE tcm_3d_data_averaging 

 821  

 822  

 823  !! 

 824  ! Description: 

 825  !  

[2680]  826  !> Define 2D output variables. 

[2353]  827  !! 

[2680]  828  SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, & 

[3241]  829  nzb_do, nzt_do ) 

[3776]  830  

[2680]  831  USE averaging, & 

 832  ONLY: diss_av, kh_av, km_av 

[2353]  833  

 834  IMPLICIT NONE 

 835  

[3083]  836  CHARACTER (LEN=*) :: grid !< name of vertical grid 

 837  CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz' 

 838  CHARACTER (LEN=*) :: variable !< name of variable 

[2353]  839  

[3129]  840  INTEGER(iwp) :: av !< flag for (non)average output 

 841  INTEGER(iwp) :: flag_nr !< number of masking flag 

 842  INTEGER(iwp) :: i !< loop index 

 843  INTEGER(iwp) :: j !< loop index 

 844  INTEGER(iwp) :: k !< loop index 

 845  INTEGER(iwp) :: nzb_do !< vertical output index (bottom) 

 846  INTEGER(iwp) :: nzt_do !< vertical output index (top) 

[2353]  847  

[3545]  848  LOGICAL :: found !< flag if output variable is found 

[3129]  849  LOGICAL :: resorted !< flag if output is already resorted 

[2353]  850  

[3545]  851  REAL(wp) :: fill_value = 9999.0_wp !< value for the _FillValue attribute 

[3004]  852  

[3014]  853  REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 

[2353]  854  !< array to which output data is resorted to 

 855  

[3129]  856  REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 

[3776]  857  

[2353]  858  found = .TRUE. 

[3129]  859  resorted = .FALSE. 

 860  ! 

 861  ! Set masking flag for topography for not resorted arrays 

 862  flag_nr = 0 

[2353]  863  

 864  SELECT CASE ( TRIM( variable ) ) 

 865  

[2680]  866  CASE ( 'diss_xy', 'diss_xz', 'diss_yz' ) 

 867  IF ( av == 0 ) THEN 

[3129]  868  to_be_resorted => diss 

[2680]  869  ELSE 

[3004]  870  IF ( .NOT. ALLOCATED( diss_av ) ) THEN 

 871  ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 872  diss_av = REAL( fill_value, KIND = wp ) 

 873  ENDIF 

[3129]  874  to_be_resorted => diss_av 

[2680]  875  ENDIF 

 876  IF ( mode == 'xy' ) grid = 'zu' 

 877  

 878  CASE ( 'kh_xy', 'kh_xz', 'kh_yz' ) 

 879  IF ( av == 0 ) THEN 

[3129]  880  to_be_resorted => kh 

[2680]  881  ELSE 

[3129]  882  IF ( .NOT. ALLOCATED( kh_av ) ) THEN 

 883  ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 884  kh_av = REAL( fill_value, KIND = wp ) 

[3004]  885  ENDIF 

[3129]  886  to_be_resorted => kh_av 

[2680]  887  ENDIF 

 888  IF ( mode == 'xy' ) grid = 'zu' 

 889  

 890  CASE ( 'km_xy', 'km_xz', 'km_yz' ) 

 891  IF ( av == 0 ) THEN 

[3129]  892  to_be_resorted => km 

[2680]  893  ELSE 

[3129]  894  IF ( .NOT. ALLOCATED( km_av ) ) THEN 

 895  ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 896  km_av = REAL( fill_value, KIND = wp ) 

[3004]  897  ENDIF 

[3129]  898  to_be_resorted => km_av 

[2680]  899  ENDIF 

 900  IF ( mode == 'xy' ) grid = 'zu' 

 901  

[2353]  902  CASE DEFAULT 

 903  found = .FALSE. 

 904  grid = 'none' 

 905  

 906  END SELECT 

[3129]  907  

 908  IF ( found .AND. .NOT. resorted ) THEN 

 909  DO i = nxl, nxr 

 910  DO j = nys, nyn 

 911  DO k = nzb_do, nzt_do 

 912  local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), & 

 913  REAL( fill_value, KIND = wp ), & 

[3776]  914  BTEST( wall_flags_0(k,j,i), flag_nr ) ) 

[3129]  915  ENDDO 

 916  ENDDO 

 917  ENDDO 

 918  ENDIF 

[3776]  919  

[2353]  920  END SUBROUTINE tcm_data_output_2d 

 921  

[3776]  922  

[2353]  923  !! 

 924  ! Description: 

 925  !  

[2680]  926  !> Define 3D output variables. 

[2353]  927  !! 

[3014]  928  SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do ) 

[2353]  929  

[3776]  930  

[2353]  931  USE averaging, & 

[2680]  932  ONLY: diss_av, kh_av, km_av 

[2353]  933  

 934  IMPLICIT NONE 

 935  

[3083]  936  CHARACTER (LEN=*) :: variable !< name of variable 

[2353]  937  

[3129]  938  INTEGER(iwp) :: av !< flag for (non)average output 

 939  INTEGER(iwp) :: flag_nr !< number of masking flag 

 940  INTEGER(iwp) :: i !< loop index 

 941  INTEGER(iwp) :: j !< loop index 

 942  INTEGER(iwp) :: k !< loop index 

 943  INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0) 

 944  INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d) 

[2353]  945  

[3129]  946  LOGICAL :: found !< flag if output variable is found 

 947  LOGICAL :: resorted !< flag if output is already resorted 

[2353]  948  

[3545]  949  REAL(wp) :: fill_value = 9999.0_wp !< value for the _FillValue attribute 

[3004]  950  

[3014]  951  REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local 

[2353]  952  !< array to which output data is resorted to 

 953  

[3129]  954  REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable 

[2353]  955  

 956  found = .TRUE. 

[3129]  957  resorted = .FALSE. 

 958  ! 

 959  ! Set masking flag for topography for not resorted arrays 

 960  flag_nr = 0 

[2353]  961  

 962  SELECT CASE ( TRIM( variable ) ) 

 963  

 964  CASE ( 'diss' ) 

 965  IF ( av == 0 ) THEN 

[3129]  966  to_be_resorted => diss 

[2353]  967  ELSE 

[3004]  968  IF ( .NOT. ALLOCATED( diss_av ) ) THEN 

 969  ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 970  diss_av = REAL( fill_value, KIND = wp ) 

 971  ENDIF 

[3129]  972  to_be_resorted => diss_av 

[2353]  973  ENDIF 

 974  

[2680]  975  CASE ( 'kh' ) 

 976  IF ( av == 0 ) THEN 

[3129]  977  to_be_resorted => kh 

[2680]  978  ELSE 

[3004]  979  IF ( .NOT. ALLOCATED( kh_av ) ) THEN 

 980  ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 981  kh_av = REAL( fill_value, KIND = wp ) 

 982  ENDIF 

[3129]  983  to_be_resorted => kh_av 

[2680]  984  ENDIF 

[2358]  985  

[2680]  986  CASE ( 'km' ) 

 987  IF ( av == 0 ) THEN 

[3129]  988  to_be_resorted => km 

[2680]  989  ELSE 

[3004]  990  IF ( .NOT. ALLOCATED( km_av ) ) THEN 

 991  ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 992  km_av = REAL( fill_value, KIND = wp ) 

 993  ENDIF 

[3129]  994  to_be_resorted => km_av 

[2680]  995  ENDIF 

[3776]  996  

[2353]  997  CASE DEFAULT 

[2680]  998  found = .FALSE. 

[2353]  999  

 1000  END SELECT 

 1001  

[3129]  1002  

 1003  IF ( found .AND. .NOT. resorted ) THEN 

 1004  DO i = nxl, nxr 

 1005  DO j = nys, nyn 

 1006  DO k = nzb_do, nzt_do 

 1007  local_pf(i,j,k) = MERGE( & 

 1008  to_be_resorted(k,j,i), & 

 1009  REAL( fill_value, KIND = wp ), & 

 1010  BTEST( wall_flags_0(k,j,i), flag_nr ) ) 

 1011  ENDDO 

 1012  ENDDO 

 1013  ENDDO 

 1014  resorted = .TRUE. 

 1015  ENDIF 

 1016  

[2680]  1017  END SUBROUTINE tcm_data_output_3d 

[2353]  1018  

 1019  

 1020  !! 

 1021  ! Description: 

 1022  !  

[2761]  1023  !> Allocate arrays and assign pointers. 

 1024  !! 

 1025  SUBROUTINE tcm_init_arrays 

 1026  

[3274]  1027  USE bulk_cloud_model_mod, & 

[2761]  1028  ONLY: collision_turbulence 

 1029  

[2938]  1030  USE pmc_interface, & 

 1031  ONLY: nested_run 

 1032  

[2761]  1033  IMPLICIT NONE 

 1034  

 1035  ALLOCATE( kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 1036  ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 1037  

 1038  ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 1039  ALLOCATE( e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 1040  ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

[3636]  1041  

[2938]  1042  ! 

[3776]  1043  ! Allocate arrays required for dissipation. 

[2938]  1044  ! Please note, if it is a nested run, arrays need to be allocated even if 

[3776]  1045  ! they do not necessarily need to be transferred, which is attributed to 

 1046  ! the design of the model coupler which allocates memory for each variable. 

[4170]  1047  ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

[3636]  1048  

[4170]  1049  IF ( rans_tke_e .OR. nested_run ) THEN 

 1050  ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 1051  ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

[2761]  1052  ENDIF 

 1053  

 1054  ! 

 1055  ! Initial assignment of pointers 

 1056  e => e_1; e_p => e_2; te_m => e_3 

 1057  

[4170]  1058  diss => diss_1 

 1059  IF ( rans_tke_e .OR. nested_run ) THEN 

[2761]  1060  diss_p => diss_2; tdiss_m => diss_3 

 1061  ENDIF 

 1062  

 1063  END SUBROUTINE tcm_init_arrays 

 1064  

 1065  

 1066  !! 

 1067  ! Description: 

 1068  !  

[2680]  1069  !> Initialization of turbulence closure module. 

[2353]  1070  !! 

 1071  SUBROUTINE tcm_init 

 1072  

 1073  USE control_parameters, & 

[3241]  1074  ONLY: bc_dirichlet_l, complex_terrain, topography 

[2353]  1075  

 1076  USE model_1d_mod, & 

[3241]  1077  ONLY: e1d, kh1d, km1d 

[2353]  1078  

 1079  IMPLICIT NONE 

 1080  

[2761]  1081  INTEGER(iwp) :: i !< loop index 

 1082  INTEGER(iwp) :: j !< loop index 

 1083  INTEGER(iwp) :: k !< loop index 

[3083]  1084  INTEGER(iwp) :: nz_s_shift !< lower shift index for scalars 

 1085  INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow 

[2353]  1086  

 1087  ! 

[2913]  1088  ! Initialize mixing length 

 1089  CALL tcm_init_mixing_length 

 1090  

 1091  ! 

[2353]  1092  ! Actions for initial runs 

 1093  IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. & 

 1094  TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN 

 1095  

 1096  IF ( INDEX( initializing_actions, 'set_1dmodel_profiles' ) /= 0 ) THEN 

[3129]  1097  

 1098  IF ( .NOT. rans_tke_e ) THEN 

[2353]  1099  ! 

[3129]  1100  ! Transfer initial profiles to the arrays of the 3D model 

 1101  DO i = nxlg, nxrg 

 1102  DO j = nysg, nyng 

 1103  e(:,j,i) = e1d 

 1104  kh(:,j,i) = kh1d 

 1105  km(:,j,i) = km1d 

 1106  ENDDO 

[2353]  1107  ENDDO 

 1108  

[3129]  1109  IF ( constant_diffusion ) THEN 

 1110  e = 0.0_wp 

 1111  ENDIF 

[2353]  1112  

[4170]  1113  diss = 0.0_wp 

 1114  

[3129]  1115  ELSE 

 1116  ! 

 1117  ! In case of TKEe closure in RANS mode, do not use e, diss, and km 

 1118  ! profiles from 1D model. Instead, initialize with constant profiles 

 1119  IF ( constant_diffusion ) THEN 

 1120  km = km_constant 

 1121  kh = km / prandtl_number 

 1122  e = 0.0_wp 

 1123  ELSEIF ( e_init > 0.0_wp ) THEN 

[2519]  1124  DO i = nxlg, nxrg 

 1125  DO j = nysg, nyng 

 1126  DO k = nzb+1, nzt 

[3129]  1127  km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init ) 

[2519]  1128  ENDDO 

 1129  ENDDO 

 1130  ENDDO 

[3129]  1131  km(nzb,:,:) = km(nzb+1,:,:) 

 1132  km(nzt+1,:,:) = km(nzt,:,:) 

 1133  kh = km / prandtl_number 

 1134  e = e_init 

 1135  ELSE 

[3294]  1136  IF ( .NOT. ocean_mode ) THEN 

[3776]  1137  kh = 0.01_wp ! there must exist an initial diffusion, because 

[3129]  1138  km = 0.01_wp ! otherwise no TKE would be produced by the 

[3545]  1139  ! production terms, as long as not yet 

 1140  ! e = (u*/cm)**2 at k=nzb+1 

[3129]  1141  ELSE 

 1142  kh = 0.00001_wp 

 1143  km = 0.00001_wp 

 1144  ENDIF 

 1145  e = 0.0_wp 

[2519]  1146  ENDIF 

[3129]  1147  

 1148  DO i = nxlg, nxrg 

 1149  DO j = nysg, nyng 

 1150  DO k = nzb+1, nzt 

 1151  diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i) 

 1152  ENDDO 

 1153  ENDDO 

 1154  ENDDO 

 1155  diss(nzb,:,:) = diss(nzb+1,:,:) 

 1156  diss(nzt+1,:,:) = diss(nzt,:,:) 

 1157  

[2353]  1158  ENDIF 

 1159  

[4170]  1160  ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. & 

[2761]  1161  INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN 

[2353]  1162  

 1163  IF ( constant_diffusion ) THEN 

[3083]  1164  km = km_constant 

 1165  kh = km / prandtl_number 

 1166  e = 0.0_wp 

[2353]  1167  ELSEIF ( e_init > 0.0_wp ) THEN 

[3083]  1168  DO i = nxlg, nxrg 

 1169  DO j = nysg, nyng 

 1170  DO k = nzb+1, nzt 

 1171  km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init ) 

 1172  ENDDO 

 1173  ENDDO 

[2353]  1174  ENDDO 

 1175  km(nzb,:,:) = km(nzb+1,:,:) 

 1176  km(nzt+1,:,:) = km(nzt,:,:) 

[3083]  1177  kh = km / prandtl_number 

 1178  e = e_init 

[2353]  1179  ELSE 

[3294]  1180  IF ( .NOT. ocean_mode ) THEN 

[3776]  1181  kh = 0.01_wp ! there must exist an initial diffusion, because 

[2353]  1182  km = 0.01_wp ! otherwise no TKE would be produced by the 

[2680]  1183  ! production terms, as long as not yet 

 1184  ! e = (u*/cm)**2 at k=nzb+1 

[2353]  1185  ELSE 

 1186  kh = 0.00001_wp 

 1187  km = 0.00001_wp 

 1188  ENDIF 

 1189  e = 0.0_wp 

 1190  ENDIF 

 1191  

[3083]  1192  IF ( rans_tke_e ) THEN 

 1193  DO i = nxlg, nxrg 

 1194  DO j = nysg, nyng 

 1195  DO k = nzb+1, nzt 

 1196  diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i) 

 1197  ENDDO 

 1198  ENDDO 

 1199  ENDDO 

 1200  diss(nzb,:,:) = diss(nzb+1,:,:) 

 1201  diss(nzt+1,:,:) = diss(nzt,:,:) 

[4170]  1202  ELSE 

 1203  diss = 0.0_wp 

[3083]  1204  ENDIF 

 1205  

[2353]  1206  ENDIF 

 1207  ! 

 1208  ! Store initial profiles for output purposes etc. 

 1209  hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 ) 

 1210  hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 ) 

 1211  ! 

 1212  ! Initialize old and new time levels. 

 1213  te_m = 0.0_wp 

 1214  e_p = e 

[2519]  1215  IF ( rans_tke_e ) THEN 

 1216  tdiss_m = 0.0_wp 

 1217  diss_p = diss 

 1218  ENDIF 

[2353]  1219  

 1220  ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. & 

 1221  TRIM( initializing_actions ) == 'cyclic_fill' ) & 

 1222  THEN 

 1223  

 1224  ! 

[2761]  1225  ! In case of complex terrain and cyclic fill method as initialization, 

[3776]  1226  ! shift initial data in the vertical direction for each point in the 

[2761]  1227  ! xyplane depending on local surface height 

 1228  IF ( complex_terrain .AND. & 

 1229  TRIM( initializing_actions ) == 'cyclic_fill' ) THEN 

 1230  DO i = nxlg, nxrg 

 1231  DO j = nysg, nyng 

[4168]  1232  nz_s_shift = topo_top_ind(j,i,0) 

[2761]  1233  

 1234  e(nz_s_shift:nzt+1,j,i) = e(0:nzt+1nz_s_shift,j,i) 

 1235  km(nz_s_shift:nzt+1,j,i) = km(0:nzt+1nz_s_shift,j,i) 

 1236  kh(nz_s_shift:nzt+1,j,i) = kh(0:nzt+1nz_s_shift,j,i) 

 1237  ENDDO 

 1238  ENDDO 

[3083]  1239  IF ( rans_tke_e ) THEN 

 1240  DO i = nxlg, nxrg 

 1241  DO j = nysg, nyng 

[4168]  1242  nz_s_shift = topo_top_ind(j,i,0) 

[3083]  1243  

 1244  diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1nz_s_shift,j,i) 

 1245  ENDDO 

 1246  ENDDO 

[4170]  1247  ELSE 

 1248  diss = 0.0_wp 

[3083]  1249  ENDIF 

[2761]  1250  ENDIF 

 1251  

 1252  ! 

[2353]  1253  ! Initialization of the turbulence recycling method 

 1254  IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 

 1255  turbulent_inflow ) THEN 

[2680]  1256  mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e 

[2353]  1257  ! 

[3776]  1258  ! In case of complex terrain, determine vertical displacement at inflow 

[2761]  1259  ! boundary and adjust mean inflow profiles 

 1260  IF ( complex_terrain ) THEN 

[3083]  1261  IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. & 

 1262  nysg <= 0 .AND. nyng >= 0 ) THEN 

[4168]  1263  nz_s_shift_l = topo_top_ind(0,0,0) 

[2761]  1264  ELSE 

 1265  nz_s_shift_l = 0 

 1266  ENDIF 

 1267  #if defined( __parallel ) 

 1268  CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, & 

 1269  MPI_MAX, comm2d, ierr) 

 1270  #else 

 1271  nz_s_shift = nz_s_shift_l 

 1272  #endif 

[3083]  1273  mean_inflow_profiles(nz_s_shift:nzt+1,5) = & 

 1274  hom_sum(0:nzt+1nz_s_shift,8,0) ! e 

[2761]  1275  ENDIF 

 1276  ! 

[2353]  1277  ! Use these mean profiles at the inflow (provided that Dirichlet 

 1278  ! conditions are used) 

[3182]  1279  IF ( bc_dirichlet_l ) THEN 

[2353]  1280  DO j = nysg, nyng 

 1281  DO k = nzb, nzt+1 

 1282  e(k,j,nxlg:1) = mean_inflow_profiles(k,5) 

 1283  ENDDO 

 1284  ENDDO 

 1285  ENDIF 

 1286  ENDIF 

 1287  ! 

 1288  ! Inside buildings set TKE back to zero 

 1289  IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. & 

 1290  topography /= 'flat' ) THEN 

 1291  ! 

[2761]  1292  ! Inside buildings set TKE back to zero. 

[3083]  1293  ! Other scalars (km, kh,...) are ignored at present, 

[2353]  1294  ! maybe revise later. 

 1295  DO i = nxlg, nxrg 

 1296  DO j = nysg, nyng 

 1297  DO k = nzb, nzt 

 1298  e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, & 

 1299  BTEST( wall_flags_0(k,j,i), 0 ) ) 

 1300  ENDDO 

 1301  ENDDO 

 1302  ENDDO 

 1303  

[3083]  1304  IF ( rans_tke_e ) THEN 

 1305  DO i = nxlg, nxrg 

 1306  DO j = nysg, nyng 

 1307  DO k = nzb, nzt 

 1308  diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp, & 

 1309  BTEST( wall_flags_0(k,j,i), 0 ) ) 

 1310  ENDDO 

 1311  ENDDO 

 1312  ENDDO 

 1313  ENDIF 

[2353]  1314  ENDIF 

 1315  ! 

 1316  ! Initialize new time levels (only done in order to set boundary values 

 1317  ! including ghost points) 

 1318  e_p = e 

 1319  ! 

 1320  ! Allthough tendency arrays are set in prognostic_equations, they have 

[3083]  1321  ! to be predefined here because there they are used (but multiplied with 0) 

 1322  ! before they are set. 

[2353]  1323  te_m = 0.0_wp 

 1324  

[3083]  1325  IF ( rans_tke_e ) THEN 

 1326  diss_p = diss 

 1327  tdiss_m = 0.0_wp 

 1328  ENDIF 

 1329  

[2353]  1330  ENDIF 

 1331  

 1332  END SUBROUTINE tcm_init 

 1333  

 1334  

[4170]  1335  !! 

[2901]  1336  ! Description: 

[4170]  1337  !  

[2901]  1338  !> Precomputation of griddependent and nearwall mixing length. 

[3299]  1339  !> @todo consider walls in horizontal direction at a distance further than a 

 1340  !> single grid point (RANS mode) 

[2353]  1341  !! 

[2901]  1342  SUBROUTINE tcm_init_mixing_length 

 1343  

 1344  USE arrays_3d, & 

[2913]  1345  ONLY: dzw, ug, vg, zu, zw 

[2901]  1346  

 1347  USE control_parameters, & 

[4170]  1348  ONLY: f, message_string, wall_adjustment, wall_adjustment_factor 

[2901]  1349  

 1350  USE grid_variables, & 

 1351  ONLY: dx, dy 

 1352  

 1353  USE indices, & 

[2905]  1354  ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, & 

 1355  nzt, wall_flags_0 

 1356  

[2901]  1357  USE kinds 

 1358  

[2916]  1359  

[2901]  1360  IMPLICIT NONE 

 1361  

[4170]  1362  INTEGER(iwp) :: dist_dx !< found distance devided by dx 

 1363  INTEGER(iwp) :: i !< index variable along x 

 1364  INTEGER(iwp) :: ii !< index variable along x 

 1365  INTEGER(iwp) :: j !< index variable along y 

 1366  INTEGER(iwp) :: k !< index variable along z 

 1367  INTEGER(iwp) :: k_max_topo !< index of maximum topography height 

 1368  INTEGER(iwp) :: kk !< index variable along z 

 1369  INTEGER(iwp) :: rad_i !< search radius in grid points along x 

 1370  INTEGER(iwp) :: rad_i_l !< possible search radius to the left 

 1371  INTEGER(iwp) :: rad_i_r !< possible search radius to the right 

 1372  INTEGER(iwp) :: rad_j !< search radius in grid points along y 

 1373  INTEGER(iwp) :: rad_j_n !< possible search radius to north 

 1374  INTEGER(iwp) :: rad_j_s !< possible search radius to south 

 1375  INTEGER(iwp) :: rad_k !< search radius in grid points along z 

 1376  INTEGER(iwp) :: rad_k_b !< search radius in grid points along negative z 

 1377  INTEGER(iwp) :: rad_k_t !< search radius in grid points along positive z 

[2901]  1378  

[2915]  1379  INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yzslice of vicinity 

 1380  

[2905]  1381  INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k) 

 1382  

[4170]  1383  REAL(wp) :: radius !< search radius in meter 

[2905]  1384  

[2901]  1385  ALLOCATE( l_grid(1:nzt) ) 

 1386  ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 

 1387  ! 

[2905]  1388  ! Initialize the mixing length in case of an LESsimulation 

 1389  IF ( .NOT. rans_mode ) THEN 

[2901]  1390  ! 

[2905]  1391  ! Compute the griddependent mixing length. 

 1392  DO k = 1, nzt 

 1393  l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp 

 1394  ENDDO 

 1395  ! 

 1396  ! Initialize nearwall mixing length l_wall only in the vertical direction 

 1397  ! for the moment, multiplication with wall_adjustment_factor further below 

 1398  l_wall(nzb,:,:) = l_grid(1) 

 1399  DO k = nzb+1, nzt 

 1400  l_wall(k,:,:) = l_grid(k) 

 1401  ENDDO 

 1402  l_wall(nzt+1,:,:) = l_grid(nzt) 

[2901]  1403  

[4170]  1404  IF ( wall_adjustment ) THEN 

 1405  

 1406  DO k = 1, nzt 

 1407  IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. & 

 1408  l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN 

 1409  WRITE( message_string, * ) 'grid anisotropy exceeds ', & 

 1410  'threshold given by only local', & 

 1411  ' &horizontal reduction of near_wall ', & 

 1412  'mixing length l_wall', & 

 1413  ' &starting from height level k = ', k, & 

 1414  '.' 

 1415  CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 ) 

 1416  EXIT 

 1417  ENDIF 

 1418  ENDDO 

[2901]  1419  ! 

[4170]  1420  ! In case of topography: limit nearwall mixing length l_wall further: 

 1421  ! Go through all points of the subdomain one by one and look for the closest 

 1422  ! surface. 

 1423  ! Is this correct in the ocean case? 

 1424  DO i = nxl, nxr 

 1425  DO j = nys, nyn 

 1426  DO k = nzb+1, nzt 

[2901]  1427  ! 

[4170]  1428  ! Check if current gridpoint belongs to the atmosphere 

 1429  IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 

[2901]  1430  ! 

[4170]  1431  ! Check for neighbouring gridpoints. 

 1432  ! Vertical distance, down 

 1433  IF ( .NOT. BTEST( wall_flags_0(k1,j,i), 0 ) ) & 

 1434  l_wall(k,j,i) = MIN( l_grid(k), zu(k)  zw(k1) ) 

[2901]  1435  ! 

[4170]  1436  ! Vertical distance, up 

 1437  IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) & 

 1438  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), zw(k)  zu(k) ) 

[2901]  1439  ! 

[4170]  1440  ! ydistance 

 1441  IF ( .NOT. BTEST( wall_flags_0(k,j1,i), 0 ) .OR. & 

 1442  .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) & 

 1443  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy ) 

[2901]  1444  ! 

[4170]  1445  ! xdistance 

 1446  IF ( .NOT. BTEST( wall_flags_0(k,j,i1), 0 ) .OR. & 

 1447  .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) & 

 1448  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx ) 

[2901]  1449  ! 

[4170]  1450  ! yzdistance (vertical edges, down) 

 1451  IF ( .NOT. BTEST( wall_flags_0(k1,j1,i), 0 ) .OR. & 

 1452  .NOT. BTEST( wall_flags_0(k1,j+1,i), 0 ) ) & 

 1453  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1454  SQRT( 0.25_wp * dy**2 + & 

 1455  ( zu(k)  zw(k1) )**2 ) ) 

[2901]  1456  ! 

[4170]  1457  ! yzdistance (vertical edges, up) 

 1458  IF ( .NOT. BTEST( wall_flags_0(k+1,j1,i), 0 ) .OR. & 

 1459  .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 ) ) & 

 1460  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1461  SQRT( 0.25_wp * dy**2 + & 

 1462  ( zw(k)  zu(k) )**2 ) ) 

[3776]  1463  ! 

[4170]  1464  ! xzdistance (vertical edges, down) 

 1465  IF ( .NOT. BTEST( wall_flags_0(k1,j,i1), 0 ) .OR. & 

 1466  .NOT. BTEST( wall_flags_0(k1,j,i+1), 0 ) ) & 

 1467  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1468  SQRT( 0.25_wp * dx**2 + & 

 1469  ( zu(k)  zw(k1) )**2 ) ) 

[2901]  1470  ! 

[4170]  1471  ! xzdistance (vertical edges, up) 

 1472  IF ( .NOT. BTEST( wall_flags_0(k+1,j,i1), 0 ) .OR. & 

 1473  .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 ) ) & 

 1474  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1475  SQRT( 0.25_wp * dx**2 + & 

 1476  ( zw(k)  zu(k) )**2 ) ) 

[2901]  1477  ! 

[4170]  1478  ! xydistance (horizontal edges) 

 1479  IF ( .NOT. BTEST( wall_flags_0(k,j1,i1), 0 ) .OR. & 

 1480  .NOT. BTEST( wall_flags_0(k,j+1,i1), 0 ) .OR. & 

 1481  .NOT. BTEST( wall_flags_0(k,j1,i+1), 0 ) .OR. & 

 1482  .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) ) & 

 1483  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1484  SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) ) 

[2901]  1485  ! 

[4170]  1486  ! xyz distance (vertical and horizontal edges, down) 

 1487  IF ( .NOT. BTEST( wall_flags_0(k1,j1,i1), 0 ) .OR. & 

 1488  .NOT. BTEST( wall_flags_0(k1,j+1,i1), 0 ) .OR. & 

 1489  .NOT. BTEST( wall_flags_0(k1,j1,i+1), 0 ) .OR. & 

 1490  .NOT. BTEST( wall_flags_0(k1,j+1,i+1), 0 ) ) & 

 1491  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1492  SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 

 1493  + ( zu(k)  zw(k1) )**2 ) ) 

[2901]  1494  ! 

[4170]  1495  ! xyz distance (vertical and horizontal edges, up) 

 1496  IF ( .NOT. BTEST( wall_flags_0(k+1,j1,i1), 0 ) .OR. & 

 1497  .NOT. BTEST( wall_flags_0(k+1,j+1,i1), 0 ) .OR. & 

 1498  .NOT. BTEST( wall_flags_0(k+1,j1,i+1), 0 ) .OR. & 

 1499  .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) ) & 

 1500  l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), & 

 1501  SQRT( 0.25_wp * ( dx**2 + dy**2 ) & 

 1502  + ( zw(k)  zu(k) )**2 ) ) 

[3776]  1503  

[4170]  1504  ENDIF 

 1505  ! 

 1506  ! Adjust mixing length by walladjustment factor and limit it by l_grid 

 1507  l_wall(k,j,i) = MIN( l_wall(k,j,i) * wall_adjustment_factor, l_grid(k) ) 

[2905]  1508  

[4170]  1509  ENDDO !k loop 

 1510  ENDDO !j loop 

 1511  ENDDO !i loop 

 1512  

 1513  ENDIF !if wall_adjustment 

 1514  

[2905]  1515  ELSE 

[2901]  1516  ! 

[4170]  1517  ! Initialize the mixing length in case of a RANS simulation 

[3083]  1518  ALLOCATE( l_black(nzb:nzt+1) ) 

[2901]  1519  

[2902]  1520  ! 

[2905]  1521  ! Calculate mixing length according to Blackadar (1962) 

[2902]  1522  IF ( f /= 0.0_wp ) THEN 

[3083]  1523  l_max = 2.7E4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / & 

 1524  ABS( f ) + 1.0E10_wp 

[2902]  1525  ELSE 

 1526  l_max = 30.0_wp 

 1527  ENDIF 

 1528  

 1529  DO k = nzb, nzt 

 1530  l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max ) 

 1531  ENDDO 

 1532  

 1533  l_black(nzt+1) = l_black(nzt) 

 1534  

[2905]  1535  ! 

[3299]  1536  ! Get height level of highest topography within local subdomain 

[4170]  1537  k_max_topo = 0 

[3299]  1538  DO i = nxlg, nxrg 

 1539  DO j = nysg, nyng 

[2910]  1540  DO k = nzb+1, nzt1 

[3299]  1541  IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND. & 

[2910]  1542  k > k_max_topo ) & 

 1543  k_max_topo = k 

 1544  ENDDO 

 1545  ENDDO 

 1546  ENDDO 

[3083]  1547  

 1548  l_wall(nzb,:,:) = l_black(nzb) 

 1549  l_wall(nzt+1,:,:) = l_black(nzt+1) 

[2910]  1550  ! 

[2905]  1551  ! Limit mixing length to either nearest wall or Blackadar mixing length. 

[3776]  1552  ! For that, analyze each grid point (i/j/k) ("analysed grid point") and 

[2905]  1553  ! search within its vicinity for the shortest distance to a wall by cal 

 1554  ! culating the distance between the analysed grid point and the "viewed 

 1555  ! grid point" if it contains a wall (belongs to topography). 

 1556  DO k = nzb+1, nzt 

[2902]  1557  

[2905]  1558  radius = l_black(k) ! radius within walls are searched 

 1559  ! 

 1560  ! Set l_wall to its default maximum value (l_back) 

 1561  l_wall(k,:,:) = radius 

 1562  

 1563  ! 

 1564  ! Compute search radius as number of grid points in all directions 

 1565  rad_i = CEILING( radius / dx ) 

 1566  rad_j = CEILING( radius / dy ) 

 1567  

 1568  DO kk = 0, nztk 

 1569  rad_k_t = kk 

 1570  ! 

 1571  ! Limit upward search radius to height of maximum topography 

[2910]  1572  IF ( zu(k+kk)zu(k) >= radius .OR. k+kk >= k_max_topo ) EXIT 

[2905]  1573  ENDDO 

 1574  

 1575  DO kk = 0, k 

 1576  rad_k_b = kk 

 1577  IF ( zu(k)zu(kkk) >= radius ) EXIT 

 1578  ENDDO 

 1579  

 1580  ! 

 1581  ! Get maximum vertical radius; necessary for defining arrays 

 1582  rad_k = MAX( rad_k_b, rad_k_t ) 

 1583  ! 

 1584  ! When analysed grid point lies above maximum topography, set search 

[3776]  1585  ! radius to 0 if the distance between the analysed grid point and max 

[2905]  1586  ! topography height is larger than the maximum search radius 

[2910]  1587  IF ( zu(krad_k_b) > zu(k_max_topo) ) rad_k_b = 0 

[2905]  1588  ! 

 1589  ! Search within vicinity only if the vertical search radius is >0 

 1590  IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 ) THEN 

 1591  

[3083]  1592  !> @note shape of vicinity is larger in z direction 

 1593  !> Shape of vicinity is two grid points larger than actual search 

 1594  !> radius in vertical direction. The first and last grid point is 

 1595  !> always set to 1 to asure correct detection of topography. See 

 1596  !> function "shortest_distance" for details. 

 1597  !> 20180316, gronemeier 

[2905]  1598  ALLOCATE( vicinity(rad_k1:rad_k+1,rad_j:rad_j,rad_i:rad_i) ) 

[2915]  1599  ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) ) 

[2905]  1600  

 1601  vicinity = 1 

 1602  

 1603  DO i = nxl, nxr 

 1604  DO j = nys, nyn 

 1605  ! 

 1606  ! Start search only if (i/j/k) belongs to atmosphere 

 1607  IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN 

 1608  ! 

 1609  ! Reset topography within vicinity 

 1610  vicinity(rad_k:rad_k,:,:) = 0 

 1611  ! 

[2909]  1612  ! Copy area surrounding analysed grid point into vicinity. 

 1613  ! First, limit size of data copied to vicinity by the domain 

 1614  ! border 

[3299]  1615  !> @note limit copied area to 1 grid point in hor. dir. 

 1616  !> Ignore walls in horizontal direction which are 

 1617  !> further away than a single grid point. This allows to 

 1618  !> only search within local subdomain without the need 

 1619  !> of global topography information. 

 1620  !> The error made by this assumption are acceptable at 

 1621  !> the moment. 

 1622  !> 20181001, gronemeier 

 1623  rad_i_l = MIN( 1, rad_i, i ) 

 1624  rad_i_r = MIN( 1, rad_i, nxi ) 

[2907]  1625  

[3299]  1626  rad_j_s = MIN( 1, rad_j, j ) 

 1627  rad_j_n = MIN( 1, rad_j, nyj ) 

[2909]  1628  

 1629  CALL copy_into_vicinity( k, j, i, & 

 1630  rad_k_b, rad_k_t, & 

 1631  rad_j_s, rad_j_n, & 

 1632  rad_i_l, rad_i_r ) 

[3299]  1633  !> @note in case of cyclic boundaries, those parts of the 

 1634  !> topography which lies beyond the domain borders but 

 1635  !> still within the search radius still needs to be 

 1636  !> copied into 'vicinity'. As the effective search 

 1637  !> radius is limited to 1 at the moment, no further 

 1638  !> copying is needed. Old implementation (prior to 

 1639  !> 20181001) had this covered but used a global array. 

 1640  !> 20181001, gronemeier 

[2907]  1641  

[2905]  1642  ! 

 1643  ! Search for walls only if there is any within vicinity 

 1644  IF ( MAXVAL( vicinity(rad_k:rad_k,:,:) ) /= 0 ) THEN 

 1645  ! 

 1646  ! Search within first half (positive x) 

 1647  dist_dx = rad_i 

 1648  DO ii = 0, dist_dx 

 1649  ! 

 1650  ! Search along vertical direction only if below 

 1651  ! maximum topography 

 1652  IF ( rad_k_t > 0 ) THEN 

 1653  ! 

 1654  ! Search for walls within octant (+++) 

[2915]  1655  vic_yz = vicinity(0:rad_k+1,0:rad_j,ii) 

[2905]  1656  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1657  shortest_distance( vic_yz, .TRUE., ii ) ) 

[2905]  1658  ! 

 1659  ! Search for walls within octant (++) 

 1660  ! Switch order of array so that the analysed grid 

 1661  ! point is always located at (0/0) (required by 

 1662  ! shortest_distance"). 

[2915]  1663  vic_yz = vicinity(0:rad_k+1,0:rad_j:1,ii) 

[2905]  1664  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1665  shortest_distance( vic_yz, .TRUE., ii ) ) 

[2905]  1666  

 1667  ENDIF 

 1668  ! 

 1669  ! Search for walls within octant (+) 

[2915]  1670  vic_yz = vicinity(0:rad_k1:1,0:rad_j:1,ii) 

[2905]  1671  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1672  shortest_distance( vic_yz, .FALSE., ii ) ) 

[2905]  1673  ! 

 1674  ! Search for walls within octant (++) 

[2915]  1675  vic_yz = vicinity(0:rad_k1:1,0:rad_j,ii) 

[2905]  1676  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1677  shortest_distance( vic_yz, .FALSE., ii ) ) 

[2905]  1678  ! 

 1679  ! Reduce search along x by already found distance 

 1680  dist_dx = CEILING( l_wall(k,j,i) / dx ) 

 1681  

 1682  ENDDO 

 1683  ! 

 1684  ! Search within second half (negative x) 

 1685  DO ii = 0, dist_dx, 1 

 1686  ! 

 1687  ! Search along vertical direction only if below 

 1688  ! maximum topography 

 1689  IF ( rad_k_t > 0 ) THEN 

 1690  ! 

 1691  ! Search for walls within octant (++) 

[2915]  1692  vic_yz = vicinity(0:rad_k+1,0:rad_j,ii) 

[2905]  1693  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1694  shortest_distance( vic_yz, .TRUE., ii ) ) 

[2905]  1695  ! 

 1696  ! Search for walls within octant (+) 

 1697  ! Switch order of array so that the analysed grid 

 1698  ! point is always located at (0/0) (required by 

 1699  ! shortest_distance"). 

[2915]  1700  vic_yz = vicinity(0:rad_k+1,0:rad_j:1,ii) 

[2905]  1701  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1702  shortest_distance( vic_yz, .TRUE., ii ) ) 

[2905]  1703  

 1704  ENDIF 

 1705  ! 

 1706  ! Search for walls within octant () 

[2915]  1707  vic_yz = vicinity(0:rad_k1:1,0:rad_j:1,ii) 

[2905]  1708  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1709  shortest_distance( vic_yz, .FALSE., ii ) ) 

[2905]  1710  ! 

 1711  ! Search for walls within octant (+) 

[2915]  1712  vic_yz = vicinity(0:rad_k1:1,0:rad_j,ii) 

[2905]  1713  l_wall(k,j,i) = MIN( l_wall(k,j,i), & 

[2915]  1714  shortest_distance( vic_yz, .FALSE., ii ) ) 

[2905]  1715  ! 

 1716  ! Reduce search along x by already found distance 

 1717  dist_dx = CEILING( l_wall(k,j,i) / dx ) 

 1718  

 1719  ENDDO 

 1720  

 1721  ENDIF !Check for any walls within vicinity 

 1722  

 1723  ELSE !Check if (i,j,k) belongs to atmosphere 

 1724  

[3083]  1725  l_wall(k,j,i) = l_black(k) 

[2905]  1726  

 1727  ENDIF 

 1728  

 1729  ENDDO !j loop 

 1730  ENDDO !i loop 

 1731  

[2911]  1732  DEALLOCATE( vicinity ) 

[2915]  1733  DEALLOCATE( vic_yz ) 

[2905]  1734  

 1735  ENDIF !check vertical size of vicinity 

 1736  

 1737  ENDDO !k loop 

 1738  

[3634]  1739  !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1)) 

 1740  

[2905]  1741  ENDIF !LES or RANS mode 

 1742  

 1743  ! 

 1744  ! Set lateral boundary conditions for l_wall 

 1745  CALL exchange_horiz( l_wall, nbgp ) 

 1746  

[3634]  1747  !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) & 

 1748  !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) 

 1749  

[2905]  1750  CONTAINS 

 1751  !! 

 1752  ! Description: 

 1753  !  

[3776]  1754  !> Calculate the shortest distance between position (i/j/k)=(0/0/0) and 

[2905]  1755  !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array' 

 1756  !> closest to the origin (0/0) of 'array'. 

 1757  !! 

[3241]  1758  REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i ) 

[2905]  1759  

 1760  IMPLICIT NONE 

 1761  

 1762  LOGICAL, INTENT(IN) :: orientation !< flag if array represents an array oriented upwards (true) or downwards (false) 

 1763  

 1764  INTEGER(iwp), INTENT(IN) :: pos_i !< x position of the yzplane 'array' 

 1765  

[3299]  1766  INTEGER(iwp) :: a !< loop index 

 1767  INTEGER(iwp) :: b !< loop index 

[2905]  1768  INTEGER(iwp) :: jj !< loop index 

 1769  

[3299]  1770  INTEGER(KIND=1) :: maximum !< maximum of array along z dimension 

 1771  

[2907]  1772  INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension 

[2905]  1773  

 1774  INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yzplane at position pos_i 

 1775  

 1776  ! 

 1777  ! Get coordinate of first maximum along vertical dimension 

[3299]  1778  ! at each y position of array (similar to function maxloc but more stable). 

 1779  DO a = 0, rad_j 

 1780  loc_k(a) = rad_k+1 

 1781  maximum = MAXVAL( array(:,a) ) 

 1782  DO b = 0, rad_k+1 

[3300]  1783  IF ( array(b,a) == maximum ) THEN 

[3299]  1784  loc_k(a) = b 

 1785  EXIT 

 1786  ENDIF 

 1787  ENDDO 

 1788  ENDDO 

[2905]  1789  ! 

 1790  ! Set distance to the default maximum value (=search radius) 

 1791  shortest_distance = radius 

 1792  ! 

 1793  ! Calculate distance between position (0/0/0) and 

 1794  ! position (pos_i/jj/loc(jj)) and only save the shortest distance. 

 1795  IF ( orientation ) THEN !if array is oriented upwards 

 1796  DO jj = 0, rad_j 

[3083]  1797  shortest_distance = & 

 1798  MIN( shortest_distance, & 

 1799  SQRT( MAX(REAL(pos_i, KIND=wp)*dx0.5_wp*dx, 0.0_wp)**2 & 

 1800  + MAX(REAL(jj, KIND=wp)*dy0.5_wp*dy, 0.0_wp)**2 & 

 1801  + MAX(zw(loc_k(jj)+k1)zu(k), 0.0_wp)**2 & 

 1802  ) & 

 1803  ) 

[2905]  1804  ENDDO 

 1805  ELSE !if array is oriented downwards 

[3083]  1806  !> @note MAX within zw required to circumvent error at domain border 

 1807  !> At the domain border, if noncyclic boundary is present, the 

 1808  !> index for zw could be 1, which will be errorneous (zw(1) does 

 1809  !> not exist). The MAX function limits the index to be at least 0. 

[2905]  1810  DO jj = 0, rad_j 

[3083]  1811  shortest_distance = & 

 1812  MIN( shortest_distance, & 

 1813  SQRT( MAX(REAL(pos_i, KIND=wp)*dx0.5_wp*dx, 0.0_wp)**2 & 

 1814  + MAX(REAL(jj, KIND=wp)*dy0.5_wp*dy, 0.0_wp)**2 & 

 1815  + MAX(zu(k)zw(MAX(kloc_k(jj),0_iwp)), 0.0_wp)**2 & 

 1816  ) & 

 1817  ) 

[2905]  1818  ENDDO 

 1819  ENDIF 

[3776]  1820  

[2905]  1821  END FUNCTION 

 1822  

[2908]  1823  !! 

 1824  ! Description: 

 1825  !  

[3776]  1826  !> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point 

[2909]  1827  !> (kp,jp,ip) containing the first bit of wall_flags_0 into the array 

 1828  !> 'vicinity'. Only copy first bit as this indicates the presence of topography. 

[2908]  1829  !! 

 1830  SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir ) 

 1831  

 1832  IMPLICIT NONE 

 1833  

 1834  INTEGER(iwp), INTENT(IN) :: il !< left loop boundary 

 1835  INTEGER(iwp), INTENT(IN) :: ip !< center position in xdirection 

 1836  INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary 

 1837  INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary 

 1838  INTEGER(iwp), INTENT(IN) :: jp !< center position in ydirection 

 1839  INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary 

 1840  INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary 

 1841  INTEGER(iwp), INTENT(IN) :: kp !< center position in zdirection 

 1842  INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary 

 1843  

 1844  INTEGER(iwp) :: i !< loop index 

 1845  INTEGER(iwp) :: j !< loop index 

 1846  INTEGER(iwp) :: k !< loop index 

 1847  

[2909]  1848  DO i = il, ir 

 1849  DO j = js, jn 

 1850  DO k = kb, kt 

[2908]  1851  vicinity(k,j,i) = MERGE( 0, 1, & 

[3299]  1852  BTEST( wall_flags_0(kp+k,jp+j,ip+i), 0 ) ) 

[2908]  1853  ENDDO 

 1854  ENDDO 

 1855  ENDDO 

 1856  

 1857  END SUBROUTINE copy_into_vicinity 

 1858  

[2901]  1859  END SUBROUTINE tcm_init_mixing_length 

 1860  

 1861  

 1862  !! 

[2353]  1863  ! Description: 

 1864  !  

[2680]  1865  !> Initialize virtual velocities used later in production_e. 

[2353]  1866  !! 

[2680]  1867  SUBROUTINE production_e_init 

[2353]  1868  

[2680]  1869  USE arrays_3d, & 

 1870  ONLY: drho_air_zw, zu 

[2353]  1871  

 1872  USE control_parameters, & 

[2680]  1873  ONLY: constant_flux_layer 

[2353]  1874  

[3145]  1875  USE surface_layer_fluxes_mod, & 

 1876  ONLY: phi_m 

 1877  

[2353]  1878  IMPLICIT NONE 

 1879  

[3120]  1880  INTEGER(iwp) :: i !< grid index xdirection 

 1881  INTEGER(iwp) :: j !< grid index ydirection 

 1882  INTEGER(iwp) :: k !< grid index zdirection 

 1883  INTEGER(iwp) :: m !< running index surface elements 

[3776]  1884  

[3145]  1885  REAL(wp) :: km_sfc !< diffusion coefficient, used to compute virtual velocities 

[2353]  1886  

[2680]  1887  IF ( constant_flux_layer ) THEN 

[2353]  1888  ! 

[2680]  1889  ! Calculate a virtual velocity at the surface in a way that the 

 1890  ! vertical velocity gradient at k = 1 (u(k+1)u_0) matches the 

 1891  ! Prandtl law (w'u'/km). This gradient is used in the TKE shear 

 1892  ! production term at k=1 (see production_e_ij). 

 1893  ! The velocity gradient has to be limited in case of too small km 

 1894  ! (otherwise the timestep may be significantly reduced by large 

 1895  ! surface winds). 

 1896  ! not available in case of noncyclic boundary conditions. 

 1897  ! Default surfaces, upwardfacing 

 1898  !$OMP PARALLEL DO PRIVATE(i,j,k,m) 

[3634]  1899  !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 

 1900  !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu) 

[2680]  1901  DO m = 1, surf_def_h(0)%ns 

[2353]  1902  

[3776]  1903  i = surf_def_h(0)%i(m) 

[2680]  1904  j = surf_def_h(0)%j(m) 

 1905  k = surf_def_h(0)%k(m) 

[2353]  1906  ! 

[3776]  1907  ! Note, calculation of u_0 and v_0 is not fully accurate, as u/v 

 1908  ! and km are not on the same grid. Actually, a further 

[2680]  1909  ! interpolation of km onto the u/vgrid is necessary. However, the 

[3120]  1910  ! effect of this error is negligible. 

[3145]  1911  km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / & 

 1912  phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) ) 

 1913  

[2680]  1914  surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * & 

[3120]  1915  drho_air_zw(k1) * & 

 1916  ( zu(k+1)  zu(k1) ) / & 

[3145]  1917  ( km_sfc + 1.0E20_wp ) 

[2680]  1918  surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * & 

[3120]  1919  drho_air_zw(k1) * & 

 1920  ( zu(k+1)  zu(k1) ) / & 

[3776]  1921  ( km_sfc + 1.0E20_wp ) 

[2353]  1922  

[2680]  1923  IF ( ABS( u(k+1,j,i)  surf_def_h(0)%u_0(m) ) > & 

 1924  ABS( u(k+1,j,i)  u(k1,j,i) ) & 

 1925  ) surf_def_h(0)%u_0(m) = u(k1,j,i) 

[2353]  1926  

[2680]  1927  IF ( ABS( v(k+1,j,i)  surf_def_h(0)%v_0(m) ) > & 

 1928  ABS( v(k+1,j,i)  v(k1,j,i) ) & 

 1929  ) surf_def_h(0)%v_0(m) = v(k1,j,i) 

 1930  

 1931  ENDDO 

[2353]  1932  ! 

[2680]  1933  ! Default surfaces, downwardfacing surfaces 

 1934  !$OMP PARALLEL DO PRIVATE(i,j,k,m) 

[3634]  1935  !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 

 1936  !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km) 

[2680]  1937  DO m = 1, surf_def_h(1)%ns 

[2353]  1938  

[3776]  1939  i = surf_def_h(1)%i(m) 

[2680]  1940  j = surf_def_h(1)%j(m) 

 1941  k = surf_def_h(1)%k(m) 

[3130]  1942  ! 

[3776]  1943  ! Note, calculation of u_0 and v_0 is not fully accurate, as u/v 

 1944  ! and km are not on the same grid. Actually, a further 

[3130]  1945  ! interpolation of km onto the u/vgrid is necessary. However, the 

 1946  ! effect of this error is negligible. 

[2680]  1947  surf_def_h(1)%u_0(m) = u(k1,j,i)  surf_def_h(1)%usws(m) * & 

 1948  drho_air_zw(k1) * & 

 1949  ( zu(k+1)  zu(k1) ) / & 

[3776]  1950  ( km(k,j,i) + 1.0E20_wp ) 

[2680]  1951  surf_def_h(1)%v_0(m) = v(k1,j,i)  surf_def_h(1)%vsws(m) * & 

 1952  drho_air_zw(k1) * & 

 1953  ( zu(k+1)  zu(k1) ) / & 

[3776]  1954  ( km(k,j,i) + 1.0E20_wp ) 

[2353]  1955  

[2680]  1956  IF ( ABS( surf_def_h(1)%u_0(m)  u(k1,j,i) ) > & 

 1957  ABS( u(k+1,j,i)  u(k1,j,i) ) & 

 1958  ) surf_def_h(1)%u_0(m) = u(k+1,j,i) 

[2353]  1959  

[2680]  1960  IF ( ABS( surf_def_h(1)%v_0(m)  v(k1,j,i) ) > & 

 1961  ABS( v(k+1,j,i)  v(k1,j,i) ) & 

 1962  ) surf_def_h(1)%v_0(m) = v(k+1,j,i) 

[2353]  1963  

[2680]  1964  ENDDO 

[2353]  1965  ! 

[2680]  1966  ! Natural surfaces, upwardfacing 

 1967  !$OMP PARALLEL DO PRIVATE(i,j,k,m) 

[3634]  1968  !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 

 1969  !$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu) 

[2680]  1970  DO m = 1, surf_lsm_h%ns 

[2353]  1971  

[3130]  1972  i = surf_lsm_h%i(m) 

[2680]  1973  j = surf_lsm_h%j(m) 

 1974  k = surf_lsm_h%k(m) 

 1975  ! 

[3776]  1976  ! Note, calculation of u_0 and v_0 is not fully accurate, as u/v 

 1977  ! and km are not on the same grid. Actually, a further 

[2680]  1978  ! interpolation of km onto the u/vgrid is necessary. However, the 

[3130]  1979  ! effect of this error is negligible. 

[3145]  1980  km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / & 

 1981  phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) ) 

 1982  

[3120]  1983  surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * & 

 1984  drho_air_zw(k1) * & 

[3145]  1985  ( zu(k+1)  zu(k1) ) / & 

[3776]  1986  ( km_sfc + 1.0E20_wp ) 

[3120]  1987  surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * & 

 1988  drho_air_zw(k1) * & 

 1989  ( zu(k+1)  zu(k1) ) / & 

[3145]  1990  ( km_sfc + 1.0E20_wp ) 

[2353]  1991  

[2680]  1992  IF ( ABS( u(k+1,j,i)  surf_lsm_h%u_0(m) ) > & 

 1993  ABS( u(k+1,j,i)  u(k1,j,i) ) & 

 1994  ) surf_lsm_h%u_0(m) = u(k1,j,i) 

 1995  

 1996  IF ( ABS( v(k+1,j,i)  surf_lsm_h%v_0(m) ) > & 

 1997  ABS( v(k+1,j,i)  v(k1,j,i) ) & 

 1998  ) surf_lsm_h%v_0(m) = v(k1,j,i) 

 1999  

 2000  ENDDO 

[2353]  2001  ! 

[2680]  2002  ! Urban surfaces, upwardfacing 

 2003  !$OMP PARALLEL DO PRIVATE(i,j,k,m) 

[3634]  2004  !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) & 

 2005  !$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu) 

[2680]  2006  DO m = 1, surf_usm_h%ns 

[2353]  2007  

[3130]  2008  i = surf_usm_h%i(m) 

[2680]  2009  j = surf_usm_h%j(m) 

 2010  k = surf_usm_h%k(m) 

[2353]  2011  ! 

[3776]  2012  ! Note, calculation of u_0 and v_0 is not fully accurate, as u/v 

 2013  ! and km are not on the same grid. Actually, a further 

[2680]  2014  ! interpolation of km onto the u/vgrid is necessary. However, the 

[3130]  2015  ! effect of this error is negligible. 

[3145]  2016  km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / & 

 2017  phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) ) 

 2018  

[3120]  2019  surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * & 

 2020  drho_air_zw(k1) * & 

 2021  ( zu(k+1)  zu(k1) ) / & 

[3145]  2022  ( km_sfc + 1.0E20_wp ) 

[3120]  2023  surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * & 

 2024  drho_air_zw(k1) * & 

 2025  ( zu(k+1)  zu(k1) ) / & 

[3145]  2026  ( km_sfc + 1.0E20_wp ) 

[2353]  2027  

[2680]  2028  IF ( ABS( u(k+1,j,i)  surf_usm_h%u_0(m) ) > & 

 2029  ABS( u(k+1,j,i)  u(k1,j,i) ) & 

 2030  ) surf_usm_h%u_0(m) = u(k1,j,i) 

[2353]  2031  

[2680]  2032  IF ( ABS( v(k+1,j,i)  surf_usm_h%v_0(m) ) > & 

 2033  ABS( v(k+1,j,i)  v(k1,j,i) ) & 

 2034  ) surf_usm_h%v_0(m) = v(k1,j,i) 

[2353]  2035  

[2519]  2036  ENDDO 

[2353]  2037  

 2038  ENDIF 

 2039  

[2680]  2040  END SUBROUTINE production_e_init 

[2353]  2041  

 2042  

[4048]  2043  !! 

 2044  ! Description: 

 2045  !  

 2046  !> Execute modulespecific actions for all grid points 

 2047  !! 

 2048  SUBROUTINE tcm_actions( location ) 

 2049  

 2050  

 2051  CHARACTER (LEN=*) :: location !< 

 2052  

 2053  ! INTEGER(iwp) :: i !< 

 2054  ! INTEGER(iwp) :: j !< 

 2055  ! INTEGER(iwp) :: k !< 

 2056  

 2057  ! 

 2058  ! Here the modulespecific actions follow 

 2059  ! No calls for single grid points are allowed at locations before and 

 2060  ! after the timestep, since these calls are not within an i,jloop 

 2061  SELECT CASE ( location ) 

 2062  

 2063  CASE ( 'before_timestep' ) 

 2064  

 2065  

 2066  CASE ( 'before_prognostic_equations' ) 

 2067  

 2068  IF ( .NOT. constant_diffusion ) CALL production_e_init 

 2069  

 2070  

 2071  CASE ( 'after_integration' ) 

 2072  

 2073  

 2074  CASE ( 'after_timestep' ) 

 2075  

 2076  

 2077  CASE ( 'utendency' ) 

 2078  

 2079  

 2080  CASE ( 'vtendency' ) 

 2081  

 2082  

 2083  CASE ( 'wtendency' ) 

 2084  

 2085  

 2086  CASE ( 'pttendency' ) 

 2087  

 2088  

 2089  CASE ( 'satendency' ) 

 2090  

 2091  

 2092  CASE ( 'etendency' ) 

 2093  

 2094  

 2095  CASE ( 'qtendency' ) 

 2096  

 2097  

 2098  CASE ( 'stendency' ) 

 2099  

 2100  

 2101  CASE DEFAULT 

 2102  CONTINUE 

 2103  

 2104  END SELECT 

 2105  

 2106  END SUBROUTINE tcm_actions 

 2107  

 2108  

 2109  !! 

 2110  ! Description: 

 2111  !  

 2112  !> Execute modulespecific actions for grid point i,j 

 2113  !! 

 2114  SUBROUTINE tcm_actions_ij( i, j, location ) 

 2115  

 2116  

 2117  CHARACTER (LEN=*) :: location 

 2118  

 2119  INTEGER(iwp) :: i 

 2120  INTEGER(iwp) :: j 

 2121  

 2122  ! 

 2123  ! Here the modulespecific actions follow 

 2124  SELECT CASE ( location ) 

 2125  

 2126  CASE ( 'utendency' ) 

 2127  

 2128  ! Next line is to avoid compiler warning about unused variables. Please remove. 

 2129  IF ( i + j < 0 ) CONTINUE 

 2130  

 2131  CASE ( 'vtendency' ) 

 2132  

 2133  

 2134  CASE ( 'wtendency' ) 

 2135  

 2136  

 2137  CASE ( 'pttendency' ) 

 2138  

 2139  

 2140  CASE ( 'satendency' ) 

 2141  

 2142  

 2143  CASE ( 'etendency' ) 

 2144  

 2145  

 2146  CASE ( 'qtendency' ) 

 2147  

 2148  

 2149  CASE ( 'stendency' ) 

 2150  

 2151  

 2152  CASE DEFAULT 

 2153  CONTINUE 

 2154  

 2155  END SELECT 

 2156  

 2157  END SUBROUTINE tcm_actions_ij 

 2158  

 2159  

[2353]  2160  !! 

 2161  ! Description: 

 2162  !  

[2680]  2163  !> Prognostic equation for subgridscale TKE and TKE dissipation rate. 

[2353]  2164  !> Vectoroptimized version 

 2165  !! 

[3386]  2166  SUBROUTINE tcm_prognostic_equations 

[2353]  2167  

 2168  USE control_parameters, & 

[3775]  2169  ONLY: scalar_advec, tsc 

[2353]  2170  

 2171  IMPLICIT NONE 

 2172  

[2680]  2173  INTEGER(iwp) :: i !< loop index 

 2174  INTEGER(iwp) :: j !< loop index 

 2175  INTEGER(iwp) :: k !< loop index 

[2353]  2176  

[2680]  2177  REAL(wp) :: sbt !< wheighting factor for subtime step 

[2353]  2178  

 2179  ! 

 2180  ! If required, compute prognostic equation for turbulent kinetic 

 2181  ! energy (TKE) 

 2182  IF ( .NOT. constant_diffusion ) THEN 

 2183  

[3724]  2184  CALL cpu_log( log_point_s(67), 'tkeequation', 'start' ) 

[2353]  2185  

 2186  sbt = tsc(2) 

 2187  IF ( .NOT. use_upstream_for_tke ) THEN 

 2188  IF ( scalar_advec == 'bcscheme' ) THEN 

 2189  

 2190  IF ( timestep_scheme(1:5) /= 'runge' ) THEN 

 2191  ! 

 2192  ! BottChlond scheme always uses Euler time step. Thus: 

 2193  sbt = 1.0_wp 

 2194  ENDIF 

 2195  tend = 0.0_wp 

 2196  CALL advec_s_bc( e, 'e' ) 

 2197  

 2198  ENDIF 

 2199  ENDIF 

 2200  

 2201  ! 

 2202  ! TKEtendency terms with no communication 

 2203  IF ( scalar_advec /= 'bcscheme' .OR. use_upstream_for_tke ) THEN 

 2204  IF ( use_upstream_for_tke ) THEN 

 2205  tend = 0.0_wp 

 2206  CALL advec_s_up( e ) 

 2207  ELSE 

[3634]  2208  !$ACC KERNELS PRESENT(tend) 

[2353]  2209  tend = 0.0_wp 

[3634]  2210  !$ACC END KERNELS 

[2353]  2211  IF ( timestep_scheme(1:5) == 'runge' ) THEN 

 2212  IF ( ws_scheme_sca ) THEN 

[4109]  2213  CALL advec_s_ws( advc_flags_s, e, 'e', & 

 2214  bc_dirichlet_l .OR. bc_radiation_l, & 

 2215  bc_dirichlet_n .OR. bc_radiation_n, & 

 2216  bc_dirichlet_r .OR. bc_radiation_r, & 

 2217  bc_dirichlet_s .OR. bc_radiation_s ) 

[2353]  2218  ELSE 

 2219  CALL advec_s_pw( e ) 

 2220  ENDIF 

 2221  ELSE 

 2222  CALL advec_s_up( e ) 

 2223  ENDIF 

 2224  ENDIF 

 2225  ENDIF 

 2226  

[3398]  2227  CALL production_e( .FALSE. ) 

[2680]  2228  

[2353]  2229  IF ( .NOT. humidity ) THEN 

[3294]  2230  IF ( ocean_mode ) THEN 

[2353]  2231  CALL diffusion_e( prho, prho_reference ) 

 2232  ELSE 

 2233  CALL diffusion_e( pt, pt_reference ) 

 2234  ENDIF 

 2235  ELSE 

 2236  CALL diffusion_e( vpt, pt_reference ) 

 2237  ENDIF 

 2238  

 2239  ! 

 2240  ! Additional sink term for flows through plant canopies 

 2241  IF ( plant_canopy ) CALL pcm_tendency( 6 ) 

 2242  

[3684]  2243  ! CALL user_actions( 'etendency' ) ToDo: find general solution for circular dependency between modules 

[2353]  2244  

 2245  ! 

 2246  ! Prognostic equation for TKE. 

 2247  ! Eliminate negative TKE values, which can occur due to numerical 

 2248  ! reasons in the course of the integration. In such cases the old TKE 

 2249  ! value is reduced by 90%. 

[3634]  2250  !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 

 2251  !$ACC PRESENT(e, tend, te_m, wall_flags_0) & 

 2252  !$ACC PRESENT(tsc(3:3)) & 

 2253  !$ACC PRESENT(e_p) 

[2353]  2254  DO i = nxl, nxr 

 2255  DO j = nys, nyn 

 2256  DO k = nzb+1, nzt 

 2257  e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 

 2258  tsc(3) * te_m(k,j,i) ) & 

 2259  ) & 

 2260  * MERGE( 1.0_wp, 0.0_wp, & 

 2261  BTEST( wall_flags_0(k,j,i), 0 ) & 

 2262  ) 

 2263  IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i) 

 2264  ENDDO 

 2265  ENDDO 

 2266  ENDDO 

 2267  

 2268  ! 

 2269  ! Calculate tendencies for the next RungeKutta step 

 2270  IF ( timestep_scheme(1:5) == 'runge' ) THEN 

 2271  IF ( intermediate_timestep_count == 1 ) THEN 

[3634]  2272  !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 

 2273  !$ACC PRESENT(tend, te_m) 

[2353]  2274  DO i = nxl, nxr 

 2275  DO j = nys, nyn 

 2276  DO k = nzb+1, nzt 

 2277  te_m(k,j,i) = tend(k,j,i) 

 2278  ENDDO 

 2279  ENDDO 

 2280  ENDDO 

 2281  ELSEIF ( intermediate_timestep_count < & 

 2282  intermediate_timestep_count_max ) THEN 

[3634]  2283  !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) & 

 2284  !$ACC PRESENT(tend, te_m) 

[2353]  2285  DO i = nxl, nxr 

 2286  DO j = nys, nyn 

 2287  DO k = nzb+1, nzt 

 2288  te_m(k,j,i) = 9.5625_wp * tend(k,j,i) & 

 2289  + 5.3125_wp * te_m(k,j,i) 

 2290  ENDDO 

 2291  ENDDO 

 2292  ENDDO 

 2293  ENDIF 

 2294  ENDIF 

 2295  

[3724]  2296  CALL cpu_log( log_point_s(67), 'tkeequation', 'stop' ) 

[2353]  2297  

[2680]  2298  ENDIF ! TKE equation 

[2353]  2299  

 2300  ! 

[2519]  2301  ! If required, compute prognostic equation for TKE dissipation rate 

[2353]  2302  IF ( rans_tke_e ) THEN 

 2303  

[3724]  2304  CALL cpu_log( log_point_s(64), 'dissequation', 'start' ) 

[2353]  2305  

 2306  sbt = tsc(2) 

 2307  IF ( .NOT. use_upstream_for_tke ) THEN 

 2308  IF ( scalar_advec == 'bcscheme' ) THEN 

 2309  

 2310  IF ( timestep_scheme(1:5) /= 'runge' ) THEN 

 2311  ! 

 2312  ! BottChlond scheme always uses Euler time step. Thus: 

 2313  sbt = 1.0_wp 

 2314  ENDIF 

 2315  tend = 0.0_wp 

 2316  CALL advec_s_bc( diss, 'diss' ) 

 2317  

 2318  ENDIF 

 2319  ENDIF 

 2320  

 2321  ! 

 2322  ! dissipationtendency terms with no communication 

 2323  IF ( scalar_advec /= 'bcscheme' .OR. use_upstream_for_tke ) THEN 

 2324  IF ( use_upstream_for_tke ) THEN 

 2325  tend = 0.0_wp 

 2326  CALL advec_s_up( diss ) 

 2327  ELSE 

 2328  tend = 0.0_wp 

 2329  IF ( timestep_scheme(1:5) == 'runge' ) THEN 

 2330  IF ( ws_scheme_sca ) THEN 

[4109]  2331  CALL advec_s_ws( advc_flags_s, diss, 'diss', & 

 2332  bc_dirichlet_l .OR. bc_radiation_l, & 

 2333  bc_dirichlet_n .OR. bc_radiation_n, & 

 2334  bc_dirichlet_r .OR. bc_radiation_r, & 

 2335  bc_dirichlet_s .OR. bc_radiation_s ) 

[2353]  2336  ELSE 

 2337  CALL advec_s_pw( diss ) 

 2338  ENDIF 

 2339  ELSE 

 2340  CALL advec_s_up( diss ) 

 2341  ENDIF 

 2342  ENDIF 

 2343  ENDIF 

[2680]  2344  ! 

 2345  ! Production of TKE dissipation rate 

[3550]  2346  CALL production_e( .TRUE. ) 

 2347  ! 

 2348  ! Diffusion term of TKE dissipation rate 

[2353]  2349  CALL diffusion_diss 

 2350  ! 

 2351  ! Additional sink term for flows through plant canopies 

[3550]  2352  ! IF ( plant_canopy ) CALL pcm_tendency( ? ) !> @todo not yet implemented 

[2353]  2353  

[3684]  2354  ! CALL user_actions( 'etendency' ) ToDo: find general solution for circular dependency between modules 

[2353]  2355  

 2356  ! 

 2357  ! Prognostic equation for TKE dissipation. 

 2358  ! Eliminate negative dissipation values, which can occur due to numerical 

 2359  ! reasons in the course of the integration. In such cases the old 

 2360  ! dissipation value is reduced by 90%. 

 2361  DO i = nxl, nxr 

 2362  DO j = nys, nyn 

 2363  DO k = nzb+1, nzt 

 2364  diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + & 

 2365  tsc(3) * tdiss_m(k,j,i) ) & 

 2366  ) & 

 2367  * MERGE( 1.0_wp, 0.0_wp, & 

 2368  BTEST( wall_flags_0(k,j,i), 0 ) & 

 2369  ) 

 2370  IF ( diss_p(k,j,i) < 0.0_wp ) & 

 2371  diss_p(k,j,i) = 0.1_wp * diss(k,j,i) 

 2372  ENDDO 

 2373  ENDDO 

 2374  ENDDO 

 2375  

 2376  ! 

 2377  ! Calculate tendencies for the next RungeKutta step 

 2378  IF ( timestep_scheme(1:5) == 'runge' ) THEN 

 2379  IF ( intermediate_timestep_count == 1 ) THEN 

 2380  DO i = nxl, nxr 

 2381  DO j = nys, nyn 

 2382  DO k = nzb+1, nzt 

 2383  tdiss_m(k,j,i) = tend(k,j,i) 

 2384  ENDDO 

 2385  ENDDO 

 2386  ENDDO 

 2387  ELSEIF ( intermediate_timestep_count < & 

 2388  intermediate_timestep_count_max ) THEN 

 2389  DO i = nxl, nxr 

 2390  DO j = nys, nyn 

 2391  DO k = nzb+1, nzt 

 2392  tdiss_m(k,j,i) = 9.5625_wp * tend(k,j,i) & 

 2393  + 5.3125_wp * tdiss_m(k,j,i) 

 2394  ENDDO 

 2395  ENDDO 

 2396  ENDDO 

 2397  ENDIF 

 2398  ENDIF 

 2399  

[3724]  2400  CALL cpu_log( log_point_s(64), 'dissequation', 'stop' ) 

[2353]  2401  

 2402  ENDIF 

 2403  

[3386]  2404  END SUBROUTINE tcm_prognostic_equations 

[2353]  2405  

 2406  

 2407  !! 

 2408  ! Description: 

 2409  !  

[2680]  2410  !> Prognostic equation for subgridscale TKE and TKE dissipation rate. 

[2353]  2411  !> Cacheoptimized version 

 2412  !! 

[3386]  2413  SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn ) 

[2353]  2414  

 2415  USE arrays_3d, & 

[3241]  2416  ONLY: diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss, & 

 2417  flux_l_e, flux_s_diss, flux_s_e 

[2353]  2418  

[2680]  2419  USE control_parameters, & 

[3241]  2420  ONLY: tsc 

[2353]  2421  

 2422  IMPLICIT NONE 

 2423  

[2358]  2424  INTEGER(iwp) :: i !< loop index x direction 

[3083]  2425  INTEGER(iwp) :: i_omp !< first loop index of iloop in prognostic_equations 

[2358]  2426  INTEGER(iwp) :: j !< loop index y direction 

 2427  INTEGER(iwp) :: k !< loop index z direction 

[3083]  2428  INTEGER(iwp) :: tn !< task number of openmp task 

[2353]  2429  

 2430  ! 

[2680]  2431  ! If required, compute prognostic equation for turbulent kinetic 

 2432  ! energy (TKE) 

 2433  IF ( .NOT. constant_diffusion ) THEN 

[2353]  2434  

 2435  ! 

[2680]  2436  ! Tendencyterms for TKE 

 2437  tend(:,j,i) = 0.0_wp 

 2438  IF ( timestep_scheme(1:5) == 'runge' & 

 2439  .AND. .NOT. use_upstream_for_tke ) THEN 

 2440  IF ( ws_scheme_sca ) THEN 

[4109]  2441  CALL advec_s_ws( advc_flags_s, & 

 2442  i, j, e, 'e', flux_s_e, diss_s_e, & 

 2443  flux_l_e, diss_l_e , i_omp, tn, & 

 2444  bc_dirichlet_l .OR. bc_radiation_l, & 

 2445  bc_dirichlet_n .OR. bc_radiation_n, & 

 2446  bc_dirichlet_r .OR. bc_radiation_r, & 

 2447  bc_dirichlet_s .OR. bc_radiation_s ) 

[2680]  2448  ELSE 

 2449  CALL advec_s_pw( i, j, e ) 

 2450  ENDIF 

 2451  ELSE 

 2452  CALL advec_s_up( i, j, e ) 

 2453  ENDIF 

[2358]  2454  

[4048]  2455  CALL production_e_ij( i, j, .FALSE. ) 

[2373]  2456  

[2680]  2457  IF ( .NOT. humidity ) THEN 

[3294]  2458  IF ( ocean_mode ) THEN 

[4048]  2459  CALL diffusion_e_ij( i, j, prho, prho_reference ) 

[2680]  2460  ELSE 

[4048]  2461  CALL diffusion_e_ij( i, j, pt, pt_reference ) 

[2680]  2462  ENDIF 

 2463  ELSE 

[4048]  2464  CALL diffusion_e_ij( i, j, vpt, pt_reference ) 

[2680]  2465  ENDIF 

[2353]  2466  

 2467  ! 

[2680]  2468  ! Additional sink term for flows through plant canopies 

 2469  IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 ) 

[2353]  2470  

[3684]  2471  ! CALL user_actions( i, j, 'etendency' ) ToDo: find general solution for circular dependency between modules 

[2353]  2472  

 2473  ! 

[2680]  2474  ! Prognostic equation for TKE. 

 2475  ! Eliminate negative TKE values, which can occur due to numerical 

 2476  ! reasons in the course of the integration. In such cases the old 

 2477  ! TKE value is reduced by 90%. 

 2478  DO k = nzb+1, nzt 

 2479  e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + & 

 2480  tsc(3) * te_m(k,j,i) ) & 

 2481  ) & 

 2482  * MERGE( 1.0_wp, 0.0_wp, & 

 2483  BTEST( wall_flags_0(k,j,i), 0 ) & 

 2484  ) 

 2485  IF ( e_p(k,j,i) <= 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i) 

 2486  ENDDO 

[2353]  2487  

 2488  ! 

[2680]  2489  ! Calculate tendencies for the next RungeKutta step 

 2490  IF ( timestep_scheme(1:5) == 'runge' ) THEN 

 2491  IF ( intermediate_timestep_count == 1 ) THEN 

 2492  DO k = nzb+1, nzt 

 2493  te_m(k,j,i) = tend(k,j,i) 

 2494  ENDDO 

 2495  ELSEIF ( intermediate_timestep_count < & 

 2496  intermediate_timestep_count_max ) THEN 

 2497  DO k = nzb+1, nzt 

 2498  te_m(k,j,i) = 9.5625_wp * tend(k,j,i) + & 

 2499  5.3125_wp * te_m(k,j,i) 

 2500  ENDDO 

 2501  ENDIF 

 2502  ENDIF 

[2353]  2503  

[2680]  2504  ENDIF ! TKE equation 

[2353]  2505  

 2506  ! 

[2680]  2507  ! If required, compute prognostic equation for TKE dissipation rate 

 2508  IF ( rans_tke_e ) THEN 

[2353]  2509  ! 

[2680]  2510  ! Tendencyterms for dissipation 

 2511  tend(:,j,i) = 0.0_wp 

 2512  IF ( timestep_scheme(1:5) == 'runge' & 

 2513  .AND. .NOT. use_upstream_for_tke ) THEN 

 2514  IF ( ws_scheme_sca ) THEN 

[4109]  2515  CALL advec_s_ws( advc_flags_s, & 

 2516  i, j, diss, 'diss', flux_s_diss, diss_s_diss, & 

 2517  flux_l_diss, diss_l_diss, i_omp, tn, & 

 2518  bc_dirichlet_l .OR. bc_radiation_l, & 

 2519  bc_dirichlet_n .OR. bc_radiation_n, & 

 2520  bc_dirichlet_r .OR. bc_radiation_r, & 

 2521  bc_dirichlet_s .OR. bc_radiation_s ) 

[2680]  2522  ELSE 

 2523  CALL advec_s_pw( i, j, diss ) 

 2524  ENDIF 

 2525  ELSE 

 2526  CALL advec_s_up( i, j, diss ) 

 2527  ENDIF 

[2358]  2528  ! 

[2680]  2529  ! Production of TKE dissipation rate 

[4048]  2530  CALL production_e_ij( i, j, .TRUE. ) 

[3083]  2531  ! 

 2532  ! Diffusion term of TKE dissipation rate 

[4048]  2533  CALL diffusion_diss_ij( i, j ) 

[2353]  2534  ! 

[2680]  2535  ! Additional sink term for flows through plant canopies 

[3550]  2536  ! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !> @todo not yet implemented 

[2353]  2537  

[3684]  2538  ! CALL user_actions( i, j, 'disstendency' ) ToDo: find general solution for circular dependency between modules 

[2353]  2539  

 2540  ! 

[2680]  2541  ! Prognostic equation for TKE dissipation 

 2542  ! Eliminate negative dissipation values, which can occur due to 

 2543  ! numerical reasons in the course of the integration. In such cases 

 2544  ! the old dissipation value is reduced by 90%. 

 2545  DO k = nzb+1, nzt 

 2546  diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + & 

 2547  tsc(3) * tdiss_m(k,j,i) ) & 

[2353]  2548  ) & 

[2680]  2549  * MERGE( 1.0_wp, 0.0_wp, & 

[2353]  2550  BTEST( wall_flags_0(k,j,i), 0 )& 

[2680]  2551  ) 

 2552  ENDDO 

[2353]  2553  

 2554  ! 

[2680]  2555  ! Calculate tendencies for the next RungeKutta step 

 2556  IF ( timestep_scheme(1:5) == 'runge' ) THEN 

 2557  IF ( intermediate_timestep_count == 1 ) THEN 

 2558  DO k = nzb+1, nzt 

 2559  tdiss_m(k,j,i) = tend(k,j,i) 

 2560  ENDDO 

 2561  ELSEIF ( intermediate_timestep_count < & 

 2562  intermediate_timestep_count_max ) THEN 

 2563  DO k = nzb+1, nzt 

 2564  tdiss_m(k,j,i) = 9.5625_wp * tend(k,j,i) + & 

 2565  5.3125_wp * tdiss_m(k,j,i) 

 2566  ENDDO 

 2567  ENDIF 

 2568  ENDIF 

[2353]  2569  

[2680]  2570  ENDIF ! dissipation equation 

[2353]  2571  

[3386]  2572  END SUBROUTINE tcm_prognostic_equations_ij 

[2353]  2573  

 2574  

 2575  !! 

 2576  ! Description: 

 2577  !  

[2680]  2578  !> Production terms (shear + buoyancy) of the TKE. 

 2579  !> Vectoroptimized version 

 2580  !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is 

 2581  !> not considered well! 

[2353]  2582  !! 

[3398]  2583  SUBROUTINE production_e( diss_production ) 

[2353]  2584  

[2680]  2585  USE arrays_3d, & 

[3274]  2586  ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner 

[2353]  2587  

[2680]  2588  USE control_parameters, & 

[3274]  2589  ONLY: cloud_droplets, constant_flux_layer, neutral, & 

[2680]  2590  rho_reference, use_single_reference_value, use_surface_fluxes, & 

 2591  use_top_fluxes 

[2353]  2592  

[2680]  2593  USE grid_variables, & 

 2594  ONLY: ddx, dx, ddy, dy 

[2353]  2595  

[3274]  2596  USE bulk_cloud_model_mod, & 

 2597  ONLY: bulk_cloud_model 

 2598  

[2680]  2599  IMPLICIT NONE 

[2353]  2600  

[3398]  2601  LOGICAL :: diss_production 

 2602  

[2680]  2603  INTEGER(iwp) :: i !< running index xdirection 

 2604  INTEGER(iwp) :: j !< running index ydirection 

 2605  INTEGER(iwp) :: k !< running index zdirection 

 2606  INTEGER(iwp) :: l !< running index for different surface type orientation 

 2607  INTEGER(iwp) :: m !< running index surface elements 

 2608  INTEGER(iwp) :: surf_e !< end index of surface elements at given ij position 

 2609  INTEGER(iwp) :: surf_s !< start index of surface elements at given ij position 

[3359]  2610  INTEGER(iwp) :: flag_nr !< number of masking flag 

[2353]  2611  

[3545]  2612  REAL(wp) :: def !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j 

[2680]  2613  REAL(wp) :: flag !< flag to mask topography 

[3545]  2614  REAL(wp) :: k1 !< temporary factor 

 2615  REAL(wp) :: k2 !< temporary factor 

[2680]  2616  REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions  used to compute shear production at surfaces 

[3545]  2617  REAL(wp) :: theta !< virtual potential temperature 

 2618  REAL(wp) :: temp !< theta * Exnerfunction 

[3776]  2619  REAL(wp) :: sign_dir !< sign of walltke flux, depending on wall orientation 

[2680]  2620  REAL(wp) :: usvs !< momentum flux u"v" 

 2621  REAL(wp) :: vsus !< momentum flux v"u" 

 2622  REAL(wp) :: wsus !< momentum flux w"u" 

 2623  REAL(wp) :: wsvs !< momentum flux w"v" 

[2353]  2624  

[3359]  2625  REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of ucomponent in xdirection 

 2626  REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of ucomponent in ydirection 

 2627  REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of ucomponent in zdirection 

 2628  REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of vcomponent in xdirection 

 2629  REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of vcomponent in ydirection 

 2630  REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of vcomponent in zdirection 

 2631  REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of wcomponent in xdirection 

 2632  REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of wcomponent in ydirection 

 2633  REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of wcomponent in zdirection 

[3398]  2634  REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary fluxarray in zdirection 

[2353]  2635  

 2636  

 2637  

 2638  ! 

[3359]  2639  ! Calculate TKE production by shear. Calculate gradients at all grid 

 2640  ! points first, gradients at surfacebounded grid points will be 

 2641  ! overwritten further below. 

[3634]  2642  !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) & 

 2643  !$ACC PRIVATE(surf_s, surf_e) & 

 2644  !$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) & 

 2645  !$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_0) & 

 2646  !$ACC PRESENT(tend) & 

 2647  !$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) & 

 2648  !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) & 

 2649  !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3)) 

[3359]  2650  DO i = nxl, nxr 

 2651  DO j = nys, nyn 

[3634]  2652  !$ACC LOOP PRIVATE(k) 

[3359]  2653  DO k = nzb+1, nzt 

[2353]  2654  

[3359]  2655  dudx(k) = ( u(k,j,i+1)  u(k,j,i) ) * ddx 

 2656  dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1)  & 

 2657  u(k,j1,i)  u(k,j1,i+1) ) * ddy 

 2658  dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1)  & 

 2659  u(k1,j,i)  u(k1,j,i+1) ) * dd2zu(k) 

[2353]  2660  

[3359]  2661  dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1)  & 

 2662  v(k,j,i1)  v(k,j+1,i1) ) * ddx 

 2663  dvdy(k) = ( v(k,j+1,i)  v(k,j,i) ) * ddy 

 2664  dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i)  & 

 2665  v(k1,j,i)  v(k1,j+1,i) ) * dd2zu(k) 

[2353]  2666  

[3359]  2667  dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k1,j,i+1)  & 

 2668  w(k,j,i1)  w(k1,j,i1) ) * ddx 

 2669  dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k1,j+1,i)  & 

 2670  w(k,j1,i)  w(k1,j1,i) ) * ddy 

 2671  dwdz(k) = ( w(k,j,i)  w(k1,j,i) ) * ddzw(k) 

 2672  

[2680]  2673  ENDDO 

[2353]  2674  

[3359]  2675  

 2676  flag_nr = 29 

 2677  

 2678  

 2679  IF ( constant_flux_layer ) THEN 

[2353]  2680  ! 

[3359]  2681  

 2682  flag_nr = 0 

 2683  

 2684  ! Position beneath wall 

 2685  ! (2)  Will allways be executed. 

 2686  ! 'bottom and wall: use u_0,v_0 and wall functions' 

[2353]  2687  ! 

[2680]  2688  ! Compute gradients at north and southfacing surfaces. 

[3359]  2689  ! First, for default surfaces, then for urban surfaces. 

[2680]  2690  ! Note, so far no natural vertical surfaces implemented 

 2691  DO l = 0, 1 

 2692  surf_s = surf_def_v(l)%start_index(j,i) 

 2693  surf_e = surf_def_v(l)%end_index(j,i) 

[3634]  2694  !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir) 

[2680]  2695  DO m = surf_s, surf_e 

 2696  k = surf_def_v(l)%k(m) 

 2697  usvs = surf_def_v(l)%mom_flux_tke(0,m) 

 2698  wsvs = surf_def_v(l)%mom_flux_tke(1,m) 

[3359]  2699  

[2680]  2700  km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 

 2701  * 0.5_wp * dy 

[2353]  2702  ! 

[2680]  2703  ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 

 2704  sign_dir = MERGE( 1.0_wp, 1.0_wp, & 

[3359]  2705  BTEST( wall_flags_0(k,j1,i), flag_nr ) ) 

 2706  dudy(k) = sign_dir * usvs / ( km_neutral + 1E10_wp ) 

 2707  dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E10_wp ) 

[2680]  2708  ENDDO 

[2353]  2709  ! 

[2680]  2710  ! Natural surfaces 

 2711  surf_s = surf_lsm_v(l)%start_index(j,i) 

 2712  surf_e = surf_lsm_v(l)%end_index(j,i) 

[3634]  2713  !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir) 

[2680]  2714  DO m = surf_s, surf_e 

 2715  k = surf_lsm_v(l)%k(m) 

 2716  usvs = surf_lsm_v(l)%mom_flux_tke(0,m) 

 2717  wsvs = surf_lsm_v(l)%mom_flux_tke(1,m) 

[3359]  2718  

[2680]  2719  km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 

 2720  * 0.5_wp * dy 

[2353]  2721  ! 

[2680]  2722  ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 

 2723  sign_dir = MERGE( 1.0_wp, 1.0_wp, & 

[3359]  2724  BTEST( wall_flags_0(k,j1,i), flag_nr ) ) 

 2725  dudy(k) = sign_dir * usvs / ( km_neutral + 1E10_wp ) 

 2726  dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E10_wp ) 

 2727  ENDDO 

[2353]  2728  ! 

[2680]  2729  ! Urban surfaces 

 2730  surf_s = surf_usm_v(l)%start_index(j,i) 

 2731  surf_e = surf_usm_v(l)%end_index(j,i) 

[3634]  2732  !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir) 

[2680]  2733  DO m = surf_s, surf_e 

 2734  k = surf_usm_v(l)%k(m) 

 2735  usvs = surf_usm_v(l)%mom_flux_tke(0,m) 

 2736  wsvs = surf_usm_v(l)%mom_flux_tke(1,m) 

[3359]  2737  

[2680]  2738  km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 

 2739  * 0.5_wp * dy 

[2353]  2740  ! 

[2680]  2741  ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 

 2742  sign_dir = MERGE( 1.0_wp, 1.0_wp, & 

[3359]  2743  BTEST( wall_flags_0(k,j1,i), flag_nr ) ) 

 2744  dudy(k) = sign_dir * usvs / ( km_neutral + 1E10_wp ) 

 2745  dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E10_wp ) 

 2746  ENDDO 

[2680]  2747  ENDDO 

[2353]  2748  ! 

[2680]  2749  ! Compute gradients at east and westfacing walls 

 2750  DO l = 2, 3 

 2751  surf_s = surf_def_v(l)%start_index(j,i) 

 2752  surf_e = surf_def_v(l)%end_index(j,i) 

[3634]  2753  !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir) 

[2680]  2754  DO m = surf_s, surf_e 

 2755  k = surf_def_v(l)%k(m) 

 2756  vsus = surf_def_v(l)%mom_flux_tke(0,m) 

 2757  wsus = surf_def_v(l)%mom_flux_tke(1,m) 

[2353]  2758  

[2680]  2759  km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 

 2760  * 0.5_wp * dx 

[2353]  2761  ! 

[2680]  2762  ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 

 2763  sign_dir = MERGE( 1.0_wp, 1.0_wp, & 

[3359]  2764  BTEST( wall_flags_0(k,j,i1), flag_nr ) ) 

 2765  dvdx(k) = sign_dir * vsus / ( km_neutral + 1E10_wp ) 

 2766  dwdx(k) = sign_dir * wsus / ( km_neutral + 1E10_wp ) 

 2767  ENDDO 

[2353]  2768  ! 

[3359]  2769  ! Natural surfaces 

[2680]  2770  surf_s = surf_lsm_v(l)%start_index(j,i) 

 2771  surf_e = surf_lsm_v(l)%end_index(j,i) 

[3634]  2772  !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir) 

[2680]  2773  DO m = surf_s, surf_e 

 2774  k = surf_lsm_v(l)%k(m) 

 2775  vsus = surf_lsm_v(l)%mom_flux_tke(0,m) 

 2776  wsus = surf_lsm_v(l)%mom_flux_tke(1,m) 

[2353]  2777  

[2680]  2778  km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 

 2779  * 0.5_wp * dx 

[2353]  2780  ! 

[2680]  2781  ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 

 2782  sign_dir = MERGE( 1.0_wp, 1.0_wp, & 

[3359]  2783  BTEST( wall_flags_0(k,j,i1), flag_nr ) ) 

 2784  dvdx(k) = sign_dir * vsus / ( km_neutral + 1E10_wp ) 

 2785  dwdx(k) = sign_dir * wsus / ( km_neutral + 1E10_wp ) 

 2786  ENDDO 

[2353]  2787  ! 

[3359]  2788  ! Urban surfaces 

[2680]  2789  surf_s = surf_usm_v(l)%start_index(j,i) 

 2790  surf_e = surf_usm_v(l)%end_index(j,i) 

[3634]  2791  !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir) 

[2680]  2792  DO m = surf_s, surf_e 

 2793  k = surf_usm_v(l)%k(m) 

 2794  vsus = surf_usm_v(l)%mom_flux_tke(0,m) 

 2795  wsus = surf_usm_v(l)%mom_flux_tke(1,m) 

[2353]  2796  

[2680]  2797  km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp & 

 2798  * 0.5_wp * dx 

[2353]  2799  ! 

[2680]  2800  ! 1.0 for rightfacing wall, 1.0 for leftfacing wall 

 2801  sign_dir = MERGE( 1.0_wp, 1.0_wp, & 

[3359]  2802  BTEST( wall_flags_0(k,j,i1), flag_nr ) ) 

 2803  dvdx(k) = sign_dir * vsus / ( km_neutral + 1E10_wp ) 

 2804  dwdx(k) = sign_dir * wsus / ( km_neutral + 1E10_wp ) 

 2805  ENDDO 

[2680]  2806  ENDDO 

[2353]  2807  ! 

[2680]  2808  ! Compute gradients at upwardfacing surfaces 

 2809  surf_s = surf_def_h(0)%start_index(j,i) 

 2810  surf_e = surf_def_h(0)%end_index(j,i) 

[3634]  2811  !$ACC LOOP PRIVATE(m, k) 

[2680]  2812  DO m = surf_s, surf_e 

 2813  k = surf_def_h(0)%k(m) 

[2353]  2814  ! 

[3359]  2815  ! Please note, actually, an interpolation of u_0 and v_0 

 2816  ! onto the grid center would be required. However, this 

[2680]  2817  ! would require several data transfers between 2Dgrid and 

[3359]  2818  ! wall type. The effect of this missing interpolation is 

[2680]  2819  ! negligible. (See also production_e_init). 

[3359]  2820  dudz(k) = ( u(k+1,j,i)  surf_def_h(0)%u_0(m) ) * dd2zu(k) 

 2821  dvdz(k) = ( v(k+1,j,i)  surf_def_h(0)%v_0(m) ) * dd2zu(k) 

 2822  

[2680]  2823  ENDDO 

[2353]  2824  ! 

[2680]  2825  ! Natural surfaces 

 2826  surf_s = surf_lsm_h%start_index(j,i) 

 2827  surf_e = surf_lsm_h%end_index(j,i) 

[3634]  2828  !$ACC LOOP PRIVATE(m, k) 

[2680]  2829  DO m = surf_s, surf_e 

 2830  k = surf_lsm_h%k(m) 

[2519]  2831  

[3359]  2832  dudz(k) = ( u(k+1,j,i)  surf_lsm_h%u_0(m) ) * dd2zu(k) 

 2833  dvdz(k) = ( v(k+1,j,i)  surf_lsm_h%v_0(m) ) * dd2zu(k) 

 2834  

[2680]  2835  ENDDO 

[2353]  2836  ! 

[2680]  2837  ! Urban surfaces 

 2838  surf_s = surf_usm_h%start_index(j,i) 

 2839  surf_e = surf_usm_h%end_index(j,i) 

[3634]  2840  !$ACC LOOP PRIVATE(m, k) 

[2680]  2841  DO m = surf_s, surf_e 

 2842  k = surf_usm_h%k(m) 

[2519]  2843  

[3359]  2844  dudz(k) = ( u(k+1,j,i)  surf_usm_h%u_0(m) ) * dd2zu(k) 

 2845  dvdz(k) = ( v(k+1,j,i)  surf_usm_h%v_0(m) ) * dd2zu(k) 

 2846  

[2680]  2847  ENDDO 

[2353]  2848  ! 

[3359]  2849  ! Compute gradients at downwardfacing walls, only for 

[2680]  2850  ! nonnatural default surfaces 

 2851  surf_s = surf_def_h(1)%start_index(j,i) 

 2852  surf_e = surf_def_h(1)%end_index(j,i) 

[3634]  2853  !$ACC LOOP PRIVATE(m, k) 

[2680]  2854  DO m = surf_s, surf_e 

 2855  k = surf_def_h(1)%k(m) 

[2519]  2856  

[3359]  2857  dudz(k) = ( surf_def_h(1)%u_0(m)  u(k1,j,i) ) * dd2zu(k) 

 2858  dvdz(k) = ( surf_def_h(1)%v_0(m)  v(k1,j,i) ) * dd2zu(k) 

[2353]  2859  

 2860  ENDDO 

 2861  

 2862  

[3359]  2863  ENDIF 

[2353]  2864  

 2865  

[3634]  2866  !$ACC LOOP PRIVATE(k, def, flag) 

[3359]  2867  DO k = nzb+1, nzt 

[2353]  2868  

[3359]  2869  def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + & 

 2870  dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 + & 

 2871  dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 + & 

 2872  2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + & 

 2873  dwdy(k)*dvdz(k) ) 

[2353]  2874  

[3359]  2875  IF ( def < 0.0_wp ) def = 0.0_wp 

[2353]  2876  

[3359]  2877  flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) ) 

[2353]  2878  

[3398]  2879  IF ( .NOT. diss_production ) THEN 

[2353]  2880  

[3550]  2881  ! Compute tendency for TKEproduction from shear 

[3398]  2882  tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag 

 2883  

 2884  ELSE 

 2885  

[3550]  2886  ! RANS mode: Compute tendency for dissipationrateproduction from shear 

[3398]  2887  tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * & 

 2888  diss(k,j,i)/( e(k,j,i) + 1.0E20_wp ) * c_1 

 2889  

 2890  ENDIF 

 2891  

[3359]  2892  ENDDO 

[2353]  2893  

 2894  

[3359]  2895  ENDDO 

 2896  ENDDO 

[2353]  2897  

 2898  ! 

[3359]  2899  ! If required, calculate TKE production by buoyancy 

 2900  IF ( .NOT. neutral ) THEN 

[2353]  2901  

[3359]  2902  IF ( .NOT. humidity ) THEN 

[2353]  2903  

[3359]  2904  IF ( ocean_mode ) THEN 

[2353]  2905  ! 

[3359]  2906  ! So far in the ocean no special treatment of density flux 

 2907  ! in the bottom and top surface layer 

 2908  DO i = nxl, nxr 

[2680]  2909  DO j = nys, nyn 

[3398]  2910  

[2680]  2911  DO k = nzb+1, nzt 

[3398]  2912  tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i)  prho(k1,j,i) ) * dd2zu(k) 

[2680]  2913  ENDDO 

[2353]  2914  ! 

[2680]  2915  ! Treatment of nearsurface grid points, at up and down 

 2916  ! ward facing surfaces 

 2917  IF ( use_surface_fluxes ) THEN 

 2918  DO l = 0, 1 

 2919  surf_s = surf_def_h(l)%start_index(j,i) 

 2920  surf_e = surf_def_h(l)%end_index(j,i) 

[2519]  2921  DO m = surf_s, surf_e 

[2680]  2922  k = surf_def_h(l)%k(m) 

[3398]  2923  tmp_flux(k) = drho_air_zw(k1) * surf_def_h(l)%shf(m) 

[2519]  2924  ENDDO 

[2680]  2925  ENDDO 

 2926  ENDIF 

[2519]  2927  

[2680]  2928  IF ( use_top_fluxes ) THEN 

 2929  surf_s = surf_def_h(2)%start_index(j,i) 

 2930  surf_e = surf_def_h(2)%end_index(j,i) 

 2931  DO m = surf_s, surf_e 

 2932  k = surf_def_h(2)%k(m) 

[3398]  2933  tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m) 

[2353]  2934  ENDDO 

[2680]  2935  ENDIF 

[2353]  2936  

[3398]  2937  IF ( .NOT. diss_production ) THEN 

 2938  

[3550]  2939  ! Compute tendency for TKEproduction from shear 

[3398]  2940  DO k = nzb+1, nzt 

 2941  flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) ) 

 2942  tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 

 2943  MERGE( rho_reference, prho(k,j,i), & 

 2944  use_single_reference_value ) ) 

 2945  ENDDO 

 2946  

 2947  ELSE 

 2948  

[3550]  2949  ! RANS mode: Compute tendency for dissipationrateproduction from shear 

[3398]  2950  DO k = nzb+1, nzt 

 2951  flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) ) 

 2952  tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 

 2953  MERGE( rho_reference, prho(k,j,i), & 

 2954  use_single_reference_value ) ) * & 

 2955  diss(k,j,i)/( e(k,j,i) + 1.0E20_wp ) * & 

 2956  c_3 

 2957  ENDDO 

 2958  

 2959  ENDIF 

 2960  

[2680]  2961  ENDDO 

[3359]  2962  ENDDO 

[2353]  2963  

[3359]  2964  ELSE ! or IF ( .NOT. ocean_mode ) THEN 

[2353]  2965  

[3634]  2966  !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 

 2967  !$ACC PRIVATE(surf_s, surf_e) & 

 2968  !$ACC PRIVATE(tmp_flux(nzb+1:nzt)) & 

 2969  !$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_0) & 

 2970  !$ACC PRESENT(tend) & 

 2971  !$ACC PRESENT(surf_def_h(0:2)) & 

 2972  !$ACC PRESENT(surf_lsm_h) & 

 2973  !$ACC PRESENT(surf_usm_h) 

[3359]  2974  DO i = nxl, nxr 

[2353]  2975  DO j = nys, nyn 

[3359]  2976  

[3634]  2977  !$ACC LOOP PRIVATE(k) 

[2353]  2978  DO k = nzb+1, nzt 

[3398]  2979  tmp_flux(k) = 1.0_wp * kh(k,j,i) * ( pt(k+1,j,i)  pt(k1,j,i) ) * dd2zu(k) 

[2353]  2980  ENDDO 

 2981  

[2680]  2982  IF ( use_surface_fluxes ) THEN 

[2353]  2983  ! 

[2680]  2984  ! Default surfaces, up and downwardfacing 

[2353]  2985  DO l = 0, 1 

 2986  surf_s = surf_def_h(l)%start_index(j,i) 

 2987  surf_e = surf_def_h(l)%end_index(j,i) 

[3634]  2988  !$ACC LOOP PRIVATE(m, k) 

[2353]  2989  DO m = surf_s, surf_e 

 2990  k = surf_def_h(l)%k(m) 

[3398]  2991  tmp_flux(k) = drho_air_zw(k1) * surf_def_h(l)%shf(m) 

[3359]  2992  ENDDO 

[2353]  2993  ENDDO 

 2994  ! 

[2680]  2995  ! Natural surfaces 

[2353]  2996  surf_s = surf_lsm_h%start_index(j,i) 

 2997  surf_e = surf_lsm_h%end_index(j,i) 

[3634]  2998  !$ACC LOOP PRIVATE(m, k) 

[2353]  2999  DO m = surf_s, surf_e 

 3000  k = surf_lsm_h%k(m) 

[3398]  3001  tmp_flux(k) = drho_air_zw(k1) * surf_lsm_h%shf(m) 

[2353]  3002  ENDDO 

 3003  ! 

[2680]  3004  ! Urban surfaces 

[2353]  3005  surf_s = surf_usm_h%start_index(j,i) 

 3006  surf_e = surf_usm_h%end_index(j,i) 

[3634]  3007  !$ACC LOOP PRIVATE(m, k) 

[2353]  3008  DO m = surf_s, surf_e 

[2680]  3009  k = surf_usm_h%k(m) 

[3398]  3010  tmp_flux(k) = drho_air_zw(k1) * surf_usm_h%shf(m) 

[3359]  3011  ENDDO 

[2680]  3012  ENDIF 

[2353]  3013  

[2680]  3014  IF ( use_top_fluxes ) THEN 

 3015  surf_s = surf_def_h(2)%start_index(j,i) 

 3016  surf_e = surf_def_h(2)%end_index(j,i) 

[3634]  3017  !$ACC LOOP PRIVATE(m, k) 

[2680]  3018  DO m = surf_s, surf_e 

 3019  k = surf_def_h(2)%k(m) 

[3398]  3020  tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m) 

[2353]  3021  ENDDO 

[2680]  3022  ENDIF 

[3359]  3023  

[3398]  3024  IF ( .NOT. diss_production ) THEN 

 3025  

[3550]  3026  ! Compute tendency for TKEproduction from shear 

[3634]  3027  !$ACC LOOP PRIVATE(k, flag) 

[3398]  3028  DO k = nzb+1, nzt 

 3029  flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) ) 

 3030  tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 

 3031  MERGE( pt_reference, pt(k,j,i), & 

 3032  use_single_reference_value ) ) 

 3033  ENDDO 

 3034  

 3035  ELSE 

 3036  

[3550]  3037  ! RANS mode: Compute tendency for dissipationrateproduction from shear 

[3398]  3038  DO k = nzb+1, nzt 

 3039  flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) ) 

 3040  tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / & 

 3041  MERGE( pt_reference, pt(k,j,i), & 

 3042  use_single_reference_value ) ) * & 

 3043  diss(k,j,i)/( e(k,j,i) + 1.0E20_wp ) * & 

 3044  c_3 

 3045  ENDDO 

 3046  

 3047  ENDIF 

 3048  

[2680]  3049  ENDDO 

[3359]  3050  ENDDO 

[2353]  3051  

[3359]  3052  ENDIF ! from IF ( .NOT. ocean_mode ) 

[2353]  3053  


