source: palm/trunk/SOURCE/turbulence_closure_mod.f90 @ 4170

Last change on this file since 4170 was 4170, checked in by gronemeier, 22 months ago

changes in turbulence_closure_mod:

  • add performance optimizations according to K. Ketelsen to diffusion_e and tcm_diffusivities_default
  • bugfix in calculating l_wall for vertical walls
  • bugfix in using l_wall in initialization (consider wall_adjustment_factor)
  • always initialize diss and save the dissipation to that array

related changes in time_integration:

  • copy diss, diss_p, tdiss_m to GPU
  • Property svn:keywords set to Id
File size: 206.4 KB
Line 
1!> @file turbulence_closure_mod.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
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
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
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!
17! Copyright 2017-2018 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: turbulence_closure_mod.f90 4170 2019-08-19 17:12:31Z gronemeier $
27! - add performance optimizations according to K. Ketelsen
28!   to diffusion_e and tcm_diffusivities_default
29! - bugfix in calculating l_wall for vertical walls
30! - bugfix in using l_wall in initialization (consider wall_adjustment_factor)
31! - always initialize diss and save the dissipation to that array
32!
33! 4168 2019-08-16 13:50:17Z suehring
34! Replace function get_topography_top_index by topo_top_ind
35!
36! 4110 2019-07-22 17:05:21Z suehring
37! pass integer flag array as well as boundary flags to WS scalar advection
38! routine
39!
40! 4109 2019-07-22 17:00:34Z suehring
41! - Modularize setting of boundary conditions for TKE and dissipation
42! - Neumann boundary condition for TKE at model top is set also in child domain
43! - Revise setting of Neumann boundary conditions at non-cyclic lateral
44!   boundaries
45! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of
46!   an implicit Dirichlet boundary condition which implied a sink of TKE
47!   at vertical walls
48!
49! 4048 2019-06-21 21:00:21Z knoop
50! write out preprocessor directives; remove tailing whitespaces
51!
52! 3775 2019-03-04 12:40:20Z gronemeier
53! removed unused variables
54!
55! 3724 2019-02-06 16:28:23Z kanani
56! Correct double-used log_point_s units
57!
58! 3719 2019-02-06 13:10:18Z kanani
59! Changed log_point to log_point_s, otherwise this overlaps with
60! 'all progn.equations' cpu measurement.
61!
62! 3684 2019-01-20 20:20:58Z knoop
63! Remove unused variable simulated_time
64!
65! 3636 2018-12-19 13:48:34Z raasch
66! nopointer option removed
67!
68! 3634 2018-12-18 12:31:28Z knoop
69! OpenACC port for SPEC
70!
71! 3550 2018-11-21 16:01:01Z gronemeier
72! - calculate diss production same in vector and cache optimization
73! - move boundary condition for e and diss to boundary_conds
74!
75! 3545 2018-11-21 11:19:41Z gronemeier
76! - Set rans_mode according to value of turbulence_closure
77! - removed debug output
78!
79! 3430 2018-10-25 13:36:23Z maronga
80! Added support for buildings in the dynamic SGS model
81!
82! 3398 2018-10-22 19:30:24Z knoop
83! Refactored production_e and production_e_ij (removed code redundancy)
84!
85! 3386 2018-10-19 16:28:22Z gronemeier
86! Renamed tcm_prognostic to tcm_prognostic_equations
87!
88! 3385 2018-10-19 14:52:29Z knoop
89! Restructured loops and ifs in production_e to ease vectorization on GPUs
90!
91! 3300 2018-10-02 14:16:54Z gronemeier
92! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall
93!   calculation
94! - removed maxloc call as this produced different results for different
95!   compiler options
96!
97! 3294 2018-10-01 02:37:10Z raasch
98! changes concerning modularization of ocean option
99!
100! 3274 2018-09-24 15:42:55Z knoop
101! Modularization of all bulk cloud physics code components
102!
103! 3245 2018-09-13 14:08:16Z knoop
104! unused variables removed, shortest_distance has wp now
105!
106! 3183 2018-07-27 14:25:55Z suehring
107! Rename variables and remove unused variable from USE statement
108!
109! 3182 2018-07-27 13:36:03Z suehring
110! Use MOST for km only in RANS mode
111!
112! 3130 2018-07-16 11:08:55Z gronemeier
113! - move boundary condition of km and kh to tcm_diffusivities
114! - calculate km at boundaries according to MOST
115! - move phi_m to surface_layer_fluxes_mod
116!
117! 3129 2018-07-16 07:45:13Z gronemeier
118! - move limitation of diss to boundary_conds
119! - move boundary conditions for e and diss to boundary_conds
120! - consider non-default surfaces in tcm_diffusivities
121! - use z_mo within surface layer instead of calculating it
122! - resort output after case select -> reduced code duplication
123! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d
124!
125! 3121 2018-07-11 18:46:49Z gronemeier
126! - created the function phi_m
127! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
128!   surfaces
129!
130! 3120 2018-07-11 18:30:57Z gronemeier
131! - changed tcm_diffusivities to tcm_diffusivities_default
132! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
133!   and tcm_diffusivities_dynamic
134!
135! 3086 2018-06-25 09:08:04Z gronemeier
136! bugfix: set rans_const_sigma(1) = 1.3
137!
138! 3083 2018-06-19 14:03:12Z gronemeier
139! - set limits of diss at the end of prognostic equations
140! - call production_e to calculate production term of diss
141! - limit change of diss to -90% to +100%
142! - remove factor 0.5 from diffusion_diss_ij
143! - rename c_m into c_0, and c_h into c_4
144! - add rans_const_c and rans_const_sigma as namelist parameters
145! - add calculation of mixing length for profile output in case of rans_tke_e
146! - changed format of annotations to comply with doxygen standards
147! - calculate and save dissipation rate during rans_tke_l mode
148! - set bc at vertical walls for e, diss, km, kh
149! - bugfix: set l_wall = 0.0 within buildings
150! - set l_wall at bottom and top boundary (rans-mode)
151! - bugfix in production term for dissipation rate
152! - bugfix in diffusion of dissipation rate
153! - disable check for 1D model if rans_tke_e is used
154! - bugfixes for initialization (rans-mode):
155!    - correction of dissipation-rate formula
156!    - calculate km based on l_wall
157!    - initialize diss if 1D model is not used
158!
159! 3045 2018-05-28 07:55:41Z Giersch
160! Error message revised
161!
162! 3014 2018-05-09 08:42:38Z maronga
163! Bugfix: nzb_do and nzt_do were not used for 3d data output
164!
165! 3004 2018-04-27 12:33:25Z Giersch
166! Further allocation checks implemented
167!
168! 2938 2018-03-27 15:52:42Z suehring
169! Further todo's
170!
171! 2936 2018-03-27 14:49:27Z gronemeier
172! - defined l_grid only within this module
173! - Moved l_wall definition from modules.f90
174! - Get level of highest topography, used to limit upward distance calculation
175! - Consider cyclic boundary conditions for mixing length calculation
176! - Moved copy of wall_flags into subarray to subroutine
177! - Implemented l_wall calculation in case of RANS simulation
178! - Moved init of l_black to tcm_init_mixing_length
179! - Moved init_mixing_length from init_grid.f90 and
180!   renamed it to tcm_init_mixing_length
181!
182! 2764 2018-01-22 09:25:36Z gronemeier
183! Bugfix: remove duplicate SAVE statements
184!
185! 2746 2018-01-15 12:06:04Z suehring
186! Move flag plant canopy to modules
187!
188! 2718 2018-01-02 08:49:38Z maronga
189! Corrected "Former revisions" section
190!
191! 2701 2017-12-15 15:40:50Z suehring
192! Changes from last commit documented
193!
194! 2698 2017-12-14 18:46:24Z suehring
195! Bugfix in get_topography_top_index
196!
197! 2696 2017-12-14 17:12:51Z kanani
198! Initial revision
199!
200!
201!
202!
203! Authors:
204! --------
205! @author Tobias Gronemeier
206! @author Hauke Wurps
207!
208! Description:
209! ------------
210!> This module contains the available turbulence closures for PALM.
211!>
212!>
213!> @todo test initialization for all possibilities
214!> @todo add OpenMP directives whereever possible
215!> @todo Check for random disturbances
216!> @note <Enter notes on the module>
217!-----------------------------------------------------------------------------!
218 MODULE turbulence_closure_mod
219
220
221    USE arrays_3d,                                                            &
222        ONLY:  diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3,   &
223               e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m,          &
224               te_m, tend, u, v, vpt, w
225
226    USE basic_constants_and_equations_mod,                                    &
227        ONLY:  g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
228
229    USE control_parameters,                                                   &
230        ONLY:  bc_dirichlet_l,                                                &
231               bc_dirichlet_n,                                                &
232               bc_dirichlet_r,                                                &
233               bc_dirichlet_s,                                                &
234               bc_radiation_l,                                                &
235               bc_radiation_n,                                                &
236               bc_radiation_r,                                                &
237               bc_radiation_s,                                                &
238               child_domain,                                                  &
239               constant_diffusion, dt_3d, e_init, humidity,                   &
240               initializing_actions, intermediate_timestep_count,             &
241               intermediate_timestep_count_max, km_constant,                  &
242               les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
243               pt_reference, rans_mode, rans_tke_e, rans_tke_l,               &
244               timestep_scheme, turbulence_closure,                           &
245               turbulent_inflow, use_upstream_for_tke, vpt_reference,         &
246               ws_scheme_sca, current_timestep_number
247
248    USE advec_ws,                                                             &
249        ONLY:  advec_s_ws
250
251    USE advec_s_bc_mod,                                                       &
252        ONLY:  advec_s_bc
253
254    USE advec_s_pw_mod,                                                       &
255        ONLY:  advec_s_pw
256
257    USE advec_s_up_mod,                                                       &
258        ONLY:  advec_s_up
259
260    USE cpulog,                                                               &
261        ONLY:  cpu_log, log_point_s
262
263    USE indices,                                                              &
264        ONLY:  advc_flags_s,                                                  &
265               nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,    &
266               topo_top_ind,                                                  &
267               wall_flags_0
268
269    USE kinds
270
271    USE ocean_mod,                                                            &
272        ONLY:  prho_reference
273
274    USE pegrid
275
276    USE plant_canopy_model_mod,                                               &
277        ONLY:  pcm_tendency
278
279    USE statistics,                                                           &
280        ONLY:  hom, hom_sum, statistic_regions
281       
282    USE surface_mod,                                                          &
283        ONLY:  bc_h,                                                          &
284               bc_v,                                                          &
285               surf_def_h,                                                    &
286               surf_def_v,                                                    &
287               surf_lsm_h,                                                    &
288               surf_lsm_v,                                                    &
289               surf_usm_h,                                                    &
290               surf_usm_v
291
292    IMPLICIT NONE
293
294
295    REAL(wp) ::  c_0                !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
296    REAL(wp) ::  c_1                !< model constant for RANS mode
297    REAL(wp) ::  c_2                !< model constant for RANS mode
298    REAL(wp) ::  c_3                !< model constant for RANS mode
299    REAL(wp) ::  c_4                !< model constant for RANS mode
300    REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
301    REAL(wp) ::  dsig_e = 1.0_wp    !< factor to calculate Ke from Km (1/sigma_e)
302    REAL(wp) ::  dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
303
304    REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
305       (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
306
307    REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
308       (/ 1.0_wp, 1.30_wp /)
309
310    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
311
312    REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid     !< geometric mean of grid sizes dx, dy, dz
313
314    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
315
316!
317!-- Public variables
318    PUBLIC c_0, rans_const_c, rans_const_sigma
319
320    SAVE
321
322    PRIVATE
323!
324!-- Public subroutines
325    PUBLIC                                                                     &
326       tcm_boundary_conds,                                                     &
327       tcm_check_parameters,                                                   &
328       tcm_check_data_output,                                                  &
329       tcm_define_netcdf_grid,                                                 &
330       tcm_init_arrays,                                                        &
331       tcm_init,                                                               &
332       tcm_actions,                                                            &
333       tcm_prognostic_equations,                                               &
334       tcm_swap_timelevel,                                                     &
335       tcm_3d_data_averaging,                                                  &
336       tcm_data_output_2d,                                                     &
337       tcm_data_output_3d,                                                     &
338       tcm_diffusivities
339
340!
341!-- PALM interfaces:
342!-- Boundary conditions for subgrid TKE and dissipation
343    INTERFACE tcm_boundary_conds
344       MODULE PROCEDURE tcm_boundary_conds
345    END INTERFACE tcm_boundary_conds
346!
347!-- Input parameter checks to be done in check_parameters
348    INTERFACE tcm_check_parameters
349       MODULE PROCEDURE tcm_check_parameters
350    END INTERFACE tcm_check_parameters
351
352!
353!-- Data output checks for 2D/3D data to be done in check_parameters
354    INTERFACE tcm_check_data_output
355       MODULE PROCEDURE tcm_check_data_output
356    END INTERFACE tcm_check_data_output
357
358!
359!-- Definition of data output quantities
360    INTERFACE tcm_define_netcdf_grid
361       MODULE PROCEDURE tcm_define_netcdf_grid
362    END INTERFACE tcm_define_netcdf_grid
363
364!
365!-- Initialization of arrays
366    INTERFACE tcm_init_arrays
367       MODULE PROCEDURE tcm_init_arrays
368    END INTERFACE tcm_init_arrays
369
370!
371!-- Initialization actions
372    INTERFACE tcm_init
373       MODULE PROCEDURE tcm_init
374    END INTERFACE tcm_init
375
376!
377!-- Location specific actions
378    INTERFACE tcm_actions
379       MODULE PROCEDURE tcm_actions
380       MODULE PROCEDURE tcm_actions_ij
381    END INTERFACE tcm_actions
382
383!
384!-- Prognostic equations for TKE and TKE dissipation rate
385    INTERFACE tcm_prognostic_equations
386       MODULE PROCEDURE tcm_prognostic_equations
387       MODULE PROCEDURE tcm_prognostic_equations_ij
388    END INTERFACE tcm_prognostic_equations
389
390!
391!-- Swapping of time levels (required for prognostic variables)
392    INTERFACE tcm_swap_timelevel
393       MODULE PROCEDURE tcm_swap_timelevel
394    END INTERFACE tcm_swap_timelevel
395
396!
397!-- Averaging of 3D data for output
398    INTERFACE tcm_3d_data_averaging
399       MODULE PROCEDURE tcm_3d_data_averaging
400    END INTERFACE tcm_3d_data_averaging
401
402!
403!-- Data output of 2D quantities
404    INTERFACE tcm_data_output_2d
405       MODULE PROCEDURE tcm_data_output_2d
406    END INTERFACE tcm_data_output_2d
407
408!
409!-- Data output of 3D data
410    INTERFACE tcm_data_output_3d
411       MODULE PROCEDURE tcm_data_output_3d
412    END INTERFACE tcm_data_output_3d
413
414!
415!-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic
416    INTERFACE tcm_diffusivities
417       MODULE PROCEDURE tcm_diffusivities
418    END INTERFACE tcm_diffusivities
419
420
421 CONTAINS
422
423!------------------------------------------------------------------------------!
424! Description:
425! ------------
426!> Check parameters routine for turbulence closure module.
427!------------------------------------------------------------------------------!
428 SUBROUTINE tcm_boundary_conds
429
430    USE pmc_interface,                                                         &
431        ONLY : rans_mode_parent
432 
433    IMPLICIT NONE
434
435    INTEGER(iwp) ::  i  !< grid index x direction
436    INTEGER(iwp) ::  j  !< grid index y direction
437    INTEGER(iwp) ::  k  !< grid index z direction
438    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
439    INTEGER(iwp) ::  m  !< running index surface elements
440!
441!-- Boundary conditions for TKE.
442    IF ( .NOT. constant_diffusion )  THEN
443!
444!--    In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls.
445!--    Note, only TKE is prognostic in this case and dissipation is only
446!--    a diagnostic quantity.
447       IF ( .NOT. rans_mode )  THEN
448!
449!--       Horizontal walls, upward- and downward-facing
450          DO  l = 0, 1
451             !$OMP PARALLEL DO PRIVATE( i, j, k )
452             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
453             !$ACC PRESENT(bc_h, e_p)
454             DO  m = 1, bc_h(l)%ns
455                i = bc_h(l)%i(m)           
456                j = bc_h(l)%j(m)
457                k = bc_h(l)%k(m)
458                e_p(k+bc_h(l)%koff,j,i) = e_p(k,j,i)
459             ENDDO
460          ENDDO
461!
462!--       Vertical walls
463          DO  l = 0, 3
464!
465!--          Note concerning missing ACC directive for this loop: Even though 
466!--          the data structure bc_v is present, it may not contain any
467!--          allocated arrays in the flat but also in a topography case,
468!--          leading to a runtime error. Therefore, omit ACC directives
469!--          for this loop, in contrast to the bc_h loop.
470             !$OMP PARALLEL DO PRIVATE( i, j, k )
471             DO  m = 1, bc_v(l)%ns
472                i = bc_v(l)%i(m)       
473                j = bc_v(l)%j(m)
474                k = bc_v(l)%k(m)
475                e_p(k,j+bc_v(l)%joff,i+bc_v(l)%ioff) = e_p(k,j,i)
476             ENDDO
477          ENDDO
478!
479!--    In RANS mode, wall function is used as boundary condition for TKE
480       ELSE
481!
482!--       Use wall function within constant-flux layer
483!--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and
484!--       are therefore not set here.
485!
486!--       Upward-facing surfaces
487!--       Default surfaces
488          DO  m = 1, surf_def_h(0)%ns
489             i = surf_def_h(0)%i(m)
490             j = surf_def_h(0)%j(m)
491             k = surf_def_h(0)%k(m)
492             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
493          ENDDO
494!
495!--       Natural surfaces
496          DO  m = 1, surf_lsm_h%ns
497             i = surf_lsm_h%i(m)
498             j = surf_lsm_h%j(m)
499             k = surf_lsm_h%k(m)
500             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
501          ENDDO
502!
503!--       Urban surfaces
504          DO  m = 1, surf_usm_h%ns
505             i = surf_usm_h%i(m)
506             j = surf_usm_h%j(m)
507             k = surf_usm_h%k(m)
508             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
509          ENDDO
510!
511!--       Vertical surfaces
512          DO  l = 0, 3
513!
514!--          Default surfaces
515             DO  m = 1, surf_def_v(l)%ns
516                i = surf_def_v(l)%i(m)
517                j = surf_def_v(l)%j(m)
518                k = surf_def_v(l)%k(m)
519                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
520             ENDDO
521!
522!--          Natural surfaces
523             DO  m = 1, surf_lsm_v(l)%ns
524                i = surf_lsm_v(l)%i(m)
525                j = surf_lsm_v(l)%j(m)
526                k = surf_lsm_v(l)%k(m)
527                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
528             ENDDO
529!
530!--          Urban surfaces
531             DO  m = 1, surf_usm_v(l)%ns
532                i = surf_usm_v(l)%i(m)
533                j = surf_usm_v(l)%j(m)
534                k = surf_usm_v(l)%k(m)
535                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
536             ENDDO
537          ENDDO
538       ENDIF
539!
540!--    Set Neumann boundary condition for TKE at model top. Do this also
541!--    in case of a nested run.
542       !$ACC KERNELS PRESENT(e_p)
543       e_p(nzt+1,:,:) = e_p(nzt,:,:)
544       !$ACC END KERNELS
545!
546!--    Nesting case: if parent operates in RANS mode and child in LES mode,
547!--    no TKE is transfered. This case, set Neumann conditions at lateral and
548!--    top child boundaries.
549!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
550!--    is treated in the nesting.
551       If ( child_domain )  THEN
552          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
553
554             e_p(nzt+1,:,:) = e_p(nzt,:,:)
555             IF ( bc_dirichlet_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
556             IF ( bc_dirichlet_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
557             IF ( bc_dirichlet_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
558             IF ( bc_dirichlet_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
559
560          ENDIF
561       ENDIF
562!
563!--    At in- and outflow boundaries also set Neumann boundary conditions
564!--    for the SGS-TKE. An exception is made for the child domain if
565!--    both parent and child operate in RANS mode. This case no
566!--    lateral Neumann boundary conditions will be set but Dirichlet
567!--    conditions will be set in the nesting.
568       IF ( .NOT. child_domain  .AND.  .NOT. rans_mode_parent  .AND.           &
569            .NOT. rans_mode )  THEN
570          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
571             e_p(:,nys-1,:) = e_p(:,nys,:)
572             IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
573          ENDIF
574          IF ( bc_dirichlet_n  .OR.  bc_radiation_n )  THEN
575             e_p(:,nyn+1,:) = e_p(:,nyn,:)
576             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 
577          ENDIF
578          IF ( bc_dirichlet_l  .OR.  bc_radiation_l )  THEN
579             e_p(:,:,nxl-1) = e_p(:,:,nxl)
580             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 
581          ENDIF
582          IF ( bc_dirichlet_r  .OR.  bc_radiation_r )  THEN
583             e_p(:,:,nxr+1) = e_p(:,:,nxr)
584             IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:) 
585          ENDIF
586       ENDIF
587    ENDIF
588
589!
590!-- Boundary conditions for TKE dissipation rate in RANS mode.
591    IF ( rans_tke_e )  THEN
592!
593!--    Use wall function within constant-flux layer
594!--    Upward-facing surfaces
595!--    Default surfaces
596       DO  m = 1, surf_def_h(0)%ns
597          i = surf_def_h(0)%i(m)
598          j = surf_def_h(0)%j(m)
599          k = surf_def_h(0)%k(m)
600          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
601                        / ( kappa * surf_def_h(0)%z_mo(m) )
602       ENDDO
603!
604!--    Natural surfaces
605       DO  m = 1, surf_lsm_h%ns
606          i = surf_lsm_h%i(m)
607          j = surf_lsm_h%j(m)
608          k = surf_lsm_h%k(m)
609          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
610                        / ( kappa * surf_lsm_h%z_mo(m) )
611       ENDDO
612!
613!--    Urban surfaces
614       DO  m = 1, surf_usm_h%ns
615          i = surf_usm_h%i(m)
616          j = surf_usm_h%j(m)
617          k = surf_usm_h%k(m)
618          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
619                        / ( kappa * surf_usm_h%z_mo(m) )
620       ENDDO
621!
622!--    Vertical surfaces
623       DO  l = 0, 3
624!
625!--       Default surfaces
626          DO  m = 1, surf_def_v(l)%ns
627             i = surf_def_v(l)%i(m)
628             j = surf_def_v(l)%j(m)
629             k = surf_def_v(l)%k(m)
630             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
631                           / ( kappa * surf_def_v(l)%z_mo(m) )
632          ENDDO
633!
634!--       Natural surfaces
635          DO  m = 1, surf_lsm_v(l)%ns
636             i = surf_lsm_v(l)%i(m)
637             j = surf_lsm_v(l)%j(m)
638             k = surf_lsm_v(l)%k(m)
639             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
640                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
641          ENDDO
642!
643!--       Urban surfaces
644          DO  m = 1, surf_usm_v(l)%ns
645             i = surf_usm_v(l)%i(m)
646             j = surf_usm_v(l)%j(m)
647             k = surf_usm_v(l)%k(m)
648             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
649                           / ( kappa * surf_usm_v(l)%z_mo(m) )
650          ENDDO
651       ENDDO
652!
653!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
654!--    minimum value
655       DO  i = nxl, nxr
656          DO  j = nys, nyn
657             DO  k = nzb, nzt+1
658                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
659                                          2.0_wp * diss(k,j,i) ), &
660                                     0.1_wp * diss(k,j,i),        &
661                                     0.0001_wp )
662             ENDDO
663          ENDDO
664       ENDDO
665
666       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
667
668    ENDIF
669
670 END SUBROUTINE tcm_boundary_conds
671 
672!------------------------------------------------------------------------------!
673! Description:
674! ------------
675!> Check parameters routine for turbulence closure module.
676!------------------------------------------------------------------------------!
677 SUBROUTINE tcm_check_parameters
678
679    USE control_parameters,                                                    &
680        ONLY:  message_string, turbulent_inflow, turbulent_outflow
681
682    IMPLICIT NONE
683
684!
685!-- Define which turbulence closure is going to be used
686    SELECT CASE ( TRIM( turbulence_closure ) )
687
688       CASE ( 'dynamic' )
689          les_dynamic = .TRUE.
690
691       CASE ( 'Moeng_Wyngaard' )
692          les_mw = .TRUE.
693
694       CASE ( 'TKE-l' )
695          rans_tke_l = .TRUE.
696          rans_mode = .TRUE.
697
698       CASE ( 'TKE-e' )
699          rans_tke_e = .TRUE.
700          rans_mode = .TRUE.
701
702       CASE DEFAULT
703          message_string = 'Unknown turbulence closure: ' //                &
704                           TRIM( turbulence_closure )
705          CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
706
707    END SELECT
708!
709!-- Set variables for RANS mode or LES mode
710    IF ( rans_mode )  THEN
711!
712!--    Assign values to constants for RANS mode
713       dsig_e    = 1.0_wp / rans_const_sigma(1)
714       dsig_diss = 1.0_wp / rans_const_sigma(2)
715
716       c_0 = rans_const_c(0)
717       c_1 = rans_const_c(1)
718       c_2 = rans_const_c(2)
719       c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
720       c_4 = rans_const_c(4)
721
722       IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
723          message_string = 'turbulent inflow/outflow is not yet '//            &
724                           'implemented for RANS mode'
725          CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
726       ENDIF
727
728       message_string = 'RANS mode is still in development! ' //               &
729                        '&Not all features of PALM are yet compatible '//      &
730                        'with RANS mode. &Use at own risk!'
731       CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 )
732
733    ELSE
734!
735!--    LES mode
736       c_0 = 0.1_wp    !according to Lilly (1967) and Deardorff (1980)
737
738       dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
739                       !of K_e which is used in RANS mode
740
741    ENDIF
742
743 END SUBROUTINE tcm_check_parameters
744
745!------------------------------------------------------------------------------!
746! Description:
747! ------------
748!> Check data output.
749!------------------------------------------------------------------------------!
750 SUBROUTINE tcm_check_data_output( var, unit )
751
752    IMPLICIT NONE
753
754    CHARACTER (LEN=*) ::  unit     !< unit of output variable
755    CHARACTER (LEN=*) ::  var      !< name of output variable
756
757
758    SELECT CASE ( TRIM( var ) )
759
760       CASE ( 'diss' )
761          unit = 'm2/s3'
762
763       CASE ( 'kh', 'km' )
764          unit = 'm2/s'
765
766       CASE DEFAULT
767          unit = 'illegal'
768
769    END SELECT
770
771 END SUBROUTINE tcm_check_data_output
772
773
774!------------------------------------------------------------------------------!
775! Description:
776! ------------
777!> Define appropriate grid for netcdf variables.
778!> It is called out from subroutine netcdf.
779!------------------------------------------------------------------------------!
780 SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
781
782    IMPLICIT NONE
783
784    CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
785    CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
786    CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
787    CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
788
789    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
790
791    found  = .TRUE.
792
793!
794!-- Check for the grid
795    SELECT CASE ( TRIM( var ) )
796
797       CASE ( 'diss', 'diss_xy', 'diss_xz', 'diss_yz' )
798          grid_x = 'x'
799          grid_y = 'y'
800          grid_z = 'zu'
801
802       CASE ( 'kh', 'kh_xy', 'kh_xz', 'kh_yz' )
803          grid_x = 'x'
804          grid_y = 'y'
805          grid_z = 'zu'
806
807       CASE ( 'km', 'km_xy', 'km_xz', 'km_yz' )
808          grid_x = 'x'
809          grid_y = 'y'
810          grid_z = 'zu'
811
812       CASE DEFAULT
813          found  = .FALSE.
814          grid_x = 'none'
815          grid_y = 'none'
816          grid_z = 'none'
817
818    END SELECT
819
820 END SUBROUTINE tcm_define_netcdf_grid
821
822
823!------------------------------------------------------------------------------!
824! Description:
825! ------------
826!> Average 3D data.
827!------------------------------------------------------------------------------!
828 SUBROUTINE tcm_3d_data_averaging( mode, variable )
829
830
831    USE averaging,                                                             &
832        ONLY:  diss_av, kh_av, km_av
833
834    USE control_parameters,                                                    &
835        ONLY:  average_count_3d
836
837    IMPLICIT NONE
838
839    CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
840    CHARACTER (LEN=*) ::  variable   !< name of variable
841
842    INTEGER(iwp) ::  i   !< loop index
843    INTEGER(iwp) ::  j   !< loop index
844    INTEGER(iwp) ::  k   !< loop index
845
846    IF ( mode == 'allocate' )  THEN
847
848       SELECT CASE ( TRIM( variable ) )
849
850          CASE ( 'diss' )
851             IF ( .NOT. ALLOCATED( diss_av ) )  THEN
852                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
853             ENDIF
854             diss_av = 0.0_wp
855
856          CASE ( 'kh' )
857             IF ( .NOT. ALLOCATED( kh_av ) )  THEN
858                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
859             ENDIF
860             kh_av = 0.0_wp
861
862          CASE ( 'km' )
863             IF ( .NOT. ALLOCATED( km_av ) )  THEN
864                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
865             ENDIF
866             km_av = 0.0_wp
867
868          CASE DEFAULT
869             CONTINUE
870
871       END SELECT
872
873    ELSEIF ( mode == 'sum' )  THEN
874
875       SELECT CASE ( TRIM( variable ) )
876
877          CASE ( 'diss' )
878             IF ( ALLOCATED( diss_av ) ) THEN
879                DO  i = nxlg, nxrg
880                   DO  j = nysg, nyng
881                      DO  k = nzb, nzt+1
882                         diss_av(k,j,i) = diss_av(k,j,i) + diss(k,j,i)
883                      ENDDO
884                   ENDDO
885                ENDDO
886             ENDIF
887
888          CASE ( 'kh' )
889             IF ( ALLOCATED( kh_av ) ) THEN
890                DO  i = nxlg, nxrg
891                   DO  j = nysg, nyng
892                      DO  k = nzb, nzt+1
893                         kh_av(k,j,i) = kh_av(k,j,i) + kh(k,j,i)
894                      ENDDO
895                   ENDDO
896                ENDDO
897             ENDIF
898
899          CASE ( 'km' )
900             IF ( ALLOCATED( km_av ) ) THEN
901                DO  i = nxlg, nxrg
902                   DO  j = nysg, nyng
903                      DO  k = nzb, nzt+1
904                         km_av(k,j,i) = km_av(k,j,i) + km(k,j,i)
905                      ENDDO
906                   ENDDO
907                ENDDO
908             ENDIF
909
910          CASE DEFAULT
911             CONTINUE
912
913       END SELECT
914
915    ELSEIF ( mode == 'average' )  THEN
916
917       SELECT CASE ( TRIM( variable ) )
918
919          CASE ( 'diss' )
920             IF ( ALLOCATED( diss_av ) ) THEN
921                DO  i = nxlg, nxrg
922                   DO  j = nysg, nyng
923                      DO  k = nzb, nzt+1
924                         diss_av(k,j,i) = diss_av(k,j,i)                       &
925                                        / REAL( average_count_3d, KIND=wp )
926                      ENDDO
927                   ENDDO
928                ENDDO
929             ENDIF
930
931          CASE ( 'kh' )
932             IF ( ALLOCATED( kh_av ) ) THEN
933                DO  i = nxlg, nxrg
934                   DO  j = nysg, nyng
935                      DO  k = nzb, nzt+1
936                         kh_av(k,j,i) = kh_av(k,j,i)                           &
937                                        / REAL( average_count_3d, KIND=wp )
938                      ENDDO
939                   ENDDO
940                ENDDO
941             ENDIF
942
943          CASE ( 'km' )
944             IF ( ALLOCATED( km_av ) ) THEN
945                DO  i = nxlg, nxrg
946                   DO  j = nysg, nyng
947                      DO  k = nzb, nzt+1
948                         km_av(k,j,i) = km_av(k,j,i)                           &
949                                        / REAL( average_count_3d, KIND=wp )
950                      ENDDO
951                   ENDDO
952                ENDDO
953             ENDIF
954
955       END SELECT
956
957    ENDIF
958
959 END SUBROUTINE tcm_3d_data_averaging
960
961
962!------------------------------------------------------------------------------!
963! Description:
964! ------------
965!> Define 2D output variables.
966!------------------------------------------------------------------------------!
967 SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
968                                nzb_do, nzt_do )
969
970    USE averaging,                                                             &
971        ONLY:  diss_av, kh_av, km_av
972
973    IMPLICIT NONE
974
975    CHARACTER (LEN=*) ::  grid       !< name of vertical grid
976    CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
977    CHARACTER (LEN=*) ::  variable   !< name of variable
978
979    INTEGER(iwp) ::  av        !< flag for (non-)average output
980    INTEGER(iwp) ::  flag_nr   !< number of masking flag
981    INTEGER(iwp) ::  i         !< loop index
982    INTEGER(iwp) ::  j         !< loop index
983    INTEGER(iwp) ::  k         !< loop index
984    INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
985    INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
986
987    LOGICAL ::  found     !< flag if output variable is found
988    LOGICAL ::  resorted  !< flag if output is already resorted
989
990    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
991
992    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
993       !< array to which output data is resorted to
994
995    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
996
997    found = .TRUE.
998    resorted = .FALSE.
999!
1000!-- Set masking flag for topography for not resorted arrays
1001    flag_nr = 0
1002
1003    SELECT CASE ( TRIM( variable ) )
1004
1005       CASE ( 'diss_xy', 'diss_xz', 'diss_yz' )
1006          IF ( av == 0 )  THEN
1007             to_be_resorted => diss
1008          ELSE
1009             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
1010                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1011                diss_av = REAL( fill_value, KIND = wp )
1012             ENDIF
1013             to_be_resorted => diss_av
1014          ENDIF
1015          IF ( mode == 'xy' ) grid = 'zu'
1016
1017       CASE ( 'kh_xy', 'kh_xz', 'kh_yz' )
1018          IF ( av == 0 )  THEN
1019             to_be_resorted => kh
1020          ELSE
1021             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
1022                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1023                kh_av = REAL( fill_value, KIND = wp )
1024             ENDIF
1025             to_be_resorted => kh_av
1026          ENDIF
1027          IF ( mode == 'xy' ) grid = 'zu'
1028
1029       CASE ( 'km_xy', 'km_xz', 'km_yz' )
1030          IF ( av == 0 )  THEN
1031             to_be_resorted => km
1032          ELSE
1033             IF ( .NOT. ALLOCATED( km_av ) ) THEN
1034                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1035                km_av = REAL( fill_value, KIND = wp )
1036             ENDIF
1037             to_be_resorted => km_av
1038          ENDIF
1039          IF ( mode == 'xy' ) grid = 'zu'
1040
1041       CASE DEFAULT
1042          found = .FALSE.
1043          grid  = 'none'
1044
1045    END SELECT
1046
1047    IF ( found .AND. .NOT. resorted )  THEN
1048       DO  i = nxl, nxr
1049          DO  j = nys, nyn
1050             DO  k = nzb_do, nzt_do
1051                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
1052                                         REAL( fill_value, KIND = wp ),        &
1053                                         BTEST( wall_flags_0(k,j,i), flag_nr ) )
1054             ENDDO
1055          ENDDO
1056       ENDDO
1057    ENDIF
1058
1059 END SUBROUTINE tcm_data_output_2d
1060
1061
1062!------------------------------------------------------------------------------!
1063! Description:
1064! ------------
1065!> Define 3D output variables.
1066!------------------------------------------------------------------------------!
1067 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
1068
1069
1070    USE averaging,                                                             &
1071        ONLY:  diss_av, kh_av, km_av
1072
1073    IMPLICIT NONE
1074
1075    CHARACTER (LEN=*) ::  variable   !< name of variable
1076
1077    INTEGER(iwp) ::  av        !< flag for (non-)average output
1078    INTEGER(iwp) ::  flag_nr   !< number of masking flag
1079    INTEGER(iwp) ::  i         !< loop index
1080    INTEGER(iwp) ::  j         !< loop index
1081    INTEGER(iwp) ::  k         !< loop index
1082    INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
1083    INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
1084
1085    LOGICAL ::  found     !< flag if output variable is found
1086    LOGICAL ::  resorted  !< flag if output is already resorted
1087
1088    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
1089
1090    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
1091       !< array to which output data is resorted to
1092
1093    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
1094
1095    found = .TRUE.
1096    resorted = .FALSE.
1097!
1098!-- Set masking flag for topography for not resorted arrays
1099    flag_nr = 0
1100
1101    SELECT CASE ( TRIM( variable ) )
1102
1103       CASE ( 'diss' )
1104          IF ( av == 0 )  THEN
1105             to_be_resorted => diss
1106          ELSE
1107             IF ( .NOT. ALLOCATED( diss_av ) ) THEN
1108                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1109                diss_av = REAL( fill_value, KIND = wp )
1110             ENDIF
1111             to_be_resorted => diss_av
1112          ENDIF
1113
1114       CASE ( 'kh' )
1115          IF ( av == 0 )  THEN
1116             to_be_resorted => kh
1117          ELSE
1118             IF ( .NOT. ALLOCATED( kh_av ) ) THEN
1119                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1120                kh_av = REAL( fill_value, KIND = wp )
1121             ENDIF
1122             to_be_resorted => kh_av
1123          ENDIF
1124
1125       CASE ( 'km' )
1126          IF ( av == 0 )  THEN
1127             to_be_resorted => km
1128          ELSE
1129             IF ( .NOT. ALLOCATED( km_av ) ) THEN
1130                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1131                km_av = REAL( fill_value, KIND = wp )
1132             ENDIF
1133             to_be_resorted => km_av
1134          ENDIF
1135
1136       CASE DEFAULT
1137          found = .FALSE.
1138
1139    END SELECT
1140
1141
1142    IF ( found .AND. .NOT. resorted )  THEN
1143       DO  i = nxl, nxr
1144          DO  j = nys, nyn
1145             DO  k = nzb_do, nzt_do
1146                local_pf(i,j,k) = MERGE(                                 &
1147                                   to_be_resorted(k,j,i),                &
1148                                   REAL( fill_value, KIND = wp ),        &
1149                                   BTEST( wall_flags_0(k,j,i), flag_nr ) )
1150             ENDDO
1151          ENDDO
1152       ENDDO
1153       resorted = .TRUE.
1154    ENDIF
1155
1156 END SUBROUTINE tcm_data_output_3d
1157
1158
1159!------------------------------------------------------------------------------!
1160! Description:
1161! ------------
1162!> Allocate arrays and assign pointers.
1163!------------------------------------------------------------------------------!
1164 SUBROUTINE tcm_init_arrays
1165
1166    USE bulk_cloud_model_mod,                                                  &
1167        ONLY:  collision_turbulence
1168
1169    USE pmc_interface,                                                         &
1170        ONLY:  nested_run
1171
1172    IMPLICIT NONE
1173
1174    ALLOCATE( kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1175    ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1176
1177    ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1178    ALLOCATE( e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1179    ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1180
1181!
1182!-- Allocate arrays required for dissipation.
1183!-- Please note, if it is a nested run, arrays need to be allocated even if
1184!-- they do not necessarily need to be transferred, which is attributed to
1185!-- the design of the model coupler which allocates memory for each variable.
1186    ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1187
1188    IF ( rans_tke_e  .OR.  nested_run )  THEN
1189       ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1190       ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1191    ENDIF
1192
1193!
1194!-- Initial assignment of pointers
1195    e  => e_1;   e_p  => e_2;   te_m  => e_3
1196
1197    diss => diss_1
1198    IF ( rans_tke_e  .OR.  nested_run )  THEN
1199       diss_p => diss_2; tdiss_m => diss_3
1200    ENDIF
1201
1202 END SUBROUTINE tcm_init_arrays
1203
1204
1205!------------------------------------------------------------------------------!
1206! Description:
1207! ------------
1208!> Initialization of turbulence closure module.
1209!------------------------------------------------------------------------------!
1210 SUBROUTINE tcm_init
1211
1212    USE control_parameters,                                                    &
1213        ONLY:  bc_dirichlet_l, complex_terrain, topography
1214
1215    USE model_1d_mod,                                                          &
1216        ONLY:  e1d, kh1d, km1d
1217
1218    IMPLICIT NONE
1219
1220    INTEGER(iwp) :: i            !< loop index
1221    INTEGER(iwp) :: j            !< loop index
1222    INTEGER(iwp) :: k            !< loop index
1223    INTEGER(iwp) :: nz_s_shift   !< lower shift index for scalars
1224    INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
1225
1226!
1227!-- Initialize mixing length
1228    CALL tcm_init_mixing_length
1229
1230!
1231!-- Actions for initial runs
1232    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1233         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1234
1235       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1236
1237          IF ( .NOT. rans_tke_e ) THEN
1238!
1239!--          Transfer initial profiles to the arrays of the 3D model
1240             DO  i = nxlg, nxrg
1241                DO  j = nysg, nyng
1242                   e(:,j,i)  = e1d
1243                   kh(:,j,i) = kh1d
1244                   km(:,j,i) = km1d
1245                ENDDO
1246             ENDDO
1247
1248             IF ( constant_diffusion )  THEN
1249                e = 0.0_wp
1250             ENDIF
1251
1252             diss = 0.0_wp
1253
1254          ELSE
1255!
1256!--          In case of TKE-e closure in RANS mode, do not use e, diss, and km
1257!--          profiles from 1D model. Instead, initialize with constant profiles
1258             IF ( constant_diffusion )  THEN
1259                km = km_constant
1260                kh = km / prandtl_number
1261                e  = 0.0_wp
1262             ELSEIF ( e_init > 0.0_wp )  THEN
1263                DO  i = nxlg, nxrg
1264                   DO  j = nysg, nyng
1265                      DO  k = nzb+1, nzt
1266                         km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1267                      ENDDO
1268                   ENDDO
1269                ENDDO
1270                km(nzb,:,:)   = km(nzb+1,:,:)
1271                km(nzt+1,:,:) = km(nzt,:,:)
1272                kh = km / prandtl_number
1273                e  = e_init
1274             ELSE
1275                IF ( .NOT. ocean_mode )  THEN
1276                   kh   = 0.01_wp   ! there must exist an initial diffusion, because
1277                   km   = 0.01_wp   ! otherwise no TKE would be produced by the
1278                                    ! production terms, as long as not yet
1279                                    ! e = (u*/cm)**2 at k=nzb+1
1280                ELSE
1281                   kh   = 0.00001_wp
1282                   km   = 0.00001_wp
1283                ENDIF
1284                e    = 0.0_wp
1285             ENDIF
1286
1287             DO  i = nxlg, nxrg
1288                DO  j = nysg, nyng
1289                   DO  k = nzb+1, nzt
1290                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
1291                   ENDDO
1292                ENDDO
1293             ENDDO
1294             diss(nzb,:,:) = diss(nzb+1,:,:)
1295             diss(nzt+1,:,:) = diss(nzt,:,:)
1296
1297          ENDIF
1298
1299       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. &
1300                INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
1301
1302          IF ( constant_diffusion )  THEN
1303             km = km_constant
1304             kh = km / prandtl_number
1305             e  = 0.0_wp
1306          ELSEIF ( e_init > 0.0_wp )  THEN
1307             DO  i = nxlg, nxrg
1308                DO  j = nysg, nyng
1309                   DO  k = nzb+1, nzt
1310                      km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1311                   ENDDO
1312                ENDDO
1313             ENDDO
1314             km(nzb,:,:)   = km(nzb+1,:,:)
1315             km(nzt+1,:,:) = km(nzt,:,:)
1316             kh = km / prandtl_number
1317             e  = e_init
1318          ELSE
1319             IF ( .NOT. ocean_mode )  THEN
1320                kh   = 0.01_wp   ! there must exist an initial diffusion, because
1321                km   = 0.01_wp   ! otherwise no TKE would be produced by the
1322                                 ! production terms, as long as not yet
1323                                 ! e = (u*/cm)**2 at k=nzb+1
1324             ELSE
1325                kh   = 0.00001_wp
1326                km   = 0.00001_wp
1327             ENDIF
1328             e    = 0.0_wp
1329          ENDIF
1330
1331          IF ( rans_tke_e )  THEN
1332             DO  i = nxlg, nxrg
1333                DO  j = nysg, nyng
1334                   DO  k = nzb+1, nzt
1335                      diss(k,j,i) = c_0**4 * e(k,j,i)**2 / km(k,j,i)
1336                   ENDDO
1337                ENDDO
1338             ENDDO
1339             diss(nzb,:,:) = diss(nzb+1,:,:)
1340             diss(nzt+1,:,:) = diss(nzt,:,:)
1341          ELSE
1342             diss = 0.0_wp
1343          ENDIF
1344
1345       ENDIF
1346!
1347!--    Store initial profiles for output purposes etc.
1348       hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
1349       hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
1350!
1351!--    Initialize old and new time levels.
1352       te_m = 0.0_wp
1353       e_p = e
1354       IF ( rans_tke_e )  THEN
1355          tdiss_m = 0.0_wp
1356          diss_p = diss
1357       ENDIF
1358
1359    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
1360             TRIM( initializing_actions ) == 'cyclic_fill' )                   &
1361    THEN
1362
1363!
1364!--    In case of complex terrain and cyclic fill method as initialization,
1365!--    shift initial data in the vertical direction for each point in the
1366!--    x-y-plane depending on local surface height
1367       IF ( complex_terrain  .AND.                                             &
1368            TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
1369          DO  i = nxlg, nxrg
1370             DO  j = nysg, nyng
1371                nz_s_shift = topo_top_ind(j,i,0)
1372
1373                e(nz_s_shift:nzt+1,j,i)  =  e(0:nzt+1-nz_s_shift,j,i)
1374                km(nz_s_shift:nzt+1,j,i) = km(0:nzt+1-nz_s_shift,j,i)
1375                kh(nz_s_shift:nzt+1,j,i) = kh(0:nzt+1-nz_s_shift,j,i)
1376             ENDDO
1377          ENDDO
1378          IF ( rans_tke_e )  THEN
1379             DO  i = nxlg, nxrg
1380                DO  j = nysg, nyng
1381                   nz_s_shift = topo_top_ind(j,i,0)
1382
1383                   diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
1384                ENDDO
1385             ENDDO
1386          ELSE
1387             diss = 0.0_wp
1388          ENDIF
1389       ENDIF
1390
1391!
1392!--    Initialization of the turbulence recycling method
1393       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
1394            turbulent_inflow )  THEN
1395          mean_inflow_profiles(:,5) = hom_sum(:,8,0)   ! e
1396!
1397!--       In case of complex terrain, determine vertical displacement at inflow
1398!--       boundary and adjust mean inflow profiles
1399          IF ( complex_terrain )  THEN
1400             IF ( nxlg <= 0 .AND. nxrg >= 0 .AND.  &
1401                  nysg <= 0 .AND. nyng >= 0        )  THEN
1402                nz_s_shift_l = topo_top_ind(0,0,0)
1403             ELSE
1404                nz_s_shift_l = 0
1405             ENDIF
1406#if defined( __parallel )
1407             CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
1408                                MPI_MAX, comm2d, ierr)
1409#else
1410             nz_s_shift = nz_s_shift_l
1411#endif
1412             mean_inflow_profiles(nz_s_shift:nzt+1,5) =  &
1413                hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
1414          ENDIF
1415!
1416!--       Use these mean profiles at the inflow (provided that Dirichlet
1417!--       conditions are used)
1418          IF ( bc_dirichlet_l )  THEN
1419             DO  j = nysg, nyng
1420                DO  k = nzb, nzt+1
1421                   e(k,j,nxlg:-1)  = mean_inflow_profiles(k,5)
1422                ENDDO
1423             ENDDO
1424          ENDIF
1425       ENDIF
1426!
1427!--    Inside buildings set TKE back to zero
1428       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
1429            topography /= 'flat' )  THEN
1430!
1431!--       Inside buildings set TKE back to zero.
1432!--       Other scalars (km, kh,...) are ignored at present,
1433!--       maybe revise later.
1434          DO  i = nxlg, nxrg
1435             DO  j = nysg, nyng
1436                DO  k = nzb, nzt
1437                   e(k,j,i)     = MERGE( e(k,j,i), 0.0_wp,                     &
1438                                         BTEST( wall_flags_0(k,j,i), 0 ) )
1439                ENDDO
1440             ENDDO
1441          ENDDO
1442
1443          IF ( rans_tke_e )  THEN
1444             DO  i = nxlg, nxrg
1445                DO  j = nysg, nyng
1446                   DO  k = nzb, nzt
1447                      diss(k,j,i)    = MERGE( diss(k,j,i), 0.0_wp,             &
1448                                              BTEST( wall_flags_0(k,j,i), 0 ) )
1449                   ENDDO
1450                ENDDO
1451             ENDDO
1452          ENDIF
1453       ENDIF
1454!
1455!--    Initialize new time levels (only done in order to set boundary values
1456!--    including ghost points)
1457       e_p = e
1458!
1459!--    Allthough tendency arrays are set in prognostic_equations, they have
1460!--    to be predefined here because there they are used (but multiplied with 0)
1461!--    before they are set.
1462       te_m = 0.0_wp
1463
1464       IF ( rans_tke_e )  THEN
1465          diss_p = diss
1466          tdiss_m = 0.0_wp
1467       ENDIF
1468
1469    ENDIF
1470
1471 END SUBROUTINE tcm_init
1472
1473
1474!------------------------------------------------------------------------------!
1475! Description:
1476! ------------
1477!> Pre-computation of grid-dependent and near-wall mixing length.
1478!> @todo consider walls in horizontal direction at a distance further than a
1479!>       single grid point (RANS mode)
1480!------------------------------------------------------------------------------!
1481 SUBROUTINE tcm_init_mixing_length
1482
1483    USE arrays_3d,                                                             &
1484        ONLY:  dzw, ug, vg, zu, zw
1485
1486    USE control_parameters,                                                    &
1487        ONLY:  f, message_string, wall_adjustment, wall_adjustment_factor
1488
1489    USE grid_variables,                                                        &
1490        ONLY:  dx, dy
1491
1492    USE indices,                                                               &
1493        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
1494               nzt, wall_flags_0
1495
1496    USE kinds
1497
1498
1499    IMPLICIT NONE
1500
1501    INTEGER(iwp) :: dist_dx     !< found distance devided by dx
1502    INTEGER(iwp) :: i           !< index variable along x
1503    INTEGER(iwp) :: ii          !< index variable along x
1504    INTEGER(iwp) :: j           !< index variable along y
1505    INTEGER(iwp) :: k           !< index variable along z
1506    INTEGER(iwp) :: k_max_topo  !< index of maximum topography height
1507    INTEGER(iwp) :: kk          !< index variable along z
1508    INTEGER(iwp) :: rad_i       !< search radius in grid points along x
1509    INTEGER(iwp) :: rad_i_l     !< possible search radius to the left
1510    INTEGER(iwp) :: rad_i_r     !< possible search radius to the right
1511    INTEGER(iwp) :: rad_j       !< search radius in grid points along y
1512    INTEGER(iwp) :: rad_j_n     !< possible search radius to north
1513    INTEGER(iwp) :: rad_j_s     !< possible search radius to south
1514    INTEGER(iwp) :: rad_k       !< search radius in grid points along z
1515    INTEGER(iwp) :: rad_k_b     !< search radius in grid points along negative z
1516    INTEGER(iwp) :: rad_k_t     !< search radius in grid points along positive z
1517
1518    INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
1519
1520    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
1521
1522    REAL(wp) :: radius          !< search radius in meter
1523
1524    ALLOCATE( l_grid(1:nzt) )
1525    ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1526!
1527!-- Initialize the mixing length in case of an LES-simulation
1528    IF ( .NOT. rans_mode )  THEN
1529!
1530!--    Compute the grid-dependent mixing length.
1531       DO  k = 1, nzt
1532          l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
1533       ENDDO
1534!
1535!--    Initialize near-wall mixing length l_wall only in the vertical direction
1536!--    for the moment, multiplication with wall_adjustment_factor further below
1537       l_wall(nzb,:,:)   = l_grid(1)
1538       DO  k = nzb+1, nzt
1539          l_wall(k,:,:)  = l_grid(k)
1540       ENDDO
1541       l_wall(nzt+1,:,:) = l_grid(nzt)
1542
1543       IF ( wall_adjustment )  THEN
1544
1545          DO  k = 1, nzt
1546             IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
1547                  l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
1548                WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
1549                                           'threshold given by only local',        &
1550                                           ' &horizontal reduction of near_wall ', &
1551                                           'mixing length l_wall',                 &
1552                                           ' &starting from height level k = ', k, &
1553                                           '.'
1554                CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
1555                EXIT
1556             ENDIF
1557          ENDDO
1558!
1559!--       In case of topography: limit near-wall mixing length l_wall further:
1560!--       Go through all points of the subdomain one by one and look for the closest
1561!--       surface.
1562!--       Is this correct in the ocean case?
1563          DO  i = nxl, nxr
1564             DO  j = nys, nyn
1565                DO  k = nzb+1, nzt
1566!
1567!--                Check if current gridpoint belongs to the atmosphere
1568                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
1569!
1570!--                   Check for neighbouring grid-points.
1571!--                   Vertical distance, down
1572                      IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )             &
1573                         l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
1574!
1575!--                   Vertical distance, up
1576                      IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )             &
1577                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), zw(k) - zu(k) )
1578!
1579!--                   y-distance
1580                      IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .OR.         &
1581                           .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )             &
1582                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
1583!
1584!--                   x-distance
1585                      IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .OR.         &
1586                           .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )             &
1587                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
1588!
1589!--                   yz-distance (vertical edges, down)
1590                      IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 )  .OR.       &
1591                           .NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 )  )          &
1592                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
1593                                              SQRT( 0.25_wp * dy**2 +            &
1594                                             ( zu(k) - zw(k-1) )**2 ) )
1595!
1596!--                   yz-distance (vertical edges, up)
1597                      IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 )  .OR.       &
1598                           .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 )  )          &
1599                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
1600                                              SQRT( 0.25_wp * dy**2 +            &
1601                                             ( zw(k) - zu(k) )**2 ) )
1602!
1603!--                   xz-distance (vertical edges, down)
1604                      IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 )  .OR.       &
1605                           .NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 )  )          &
1606                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
1607                                              SQRT( 0.25_wp * dx**2 +            &
1608                                             ( zu(k) - zw(k-1) )**2 ) )
1609!
1610!--                   xz-distance (vertical edges, up)
1611                      IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 )  .OR.       &
1612                           .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 )  )          &
1613                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
1614                                              SQRT( 0.25_wp * dx**2 +            &
1615                                             ( zw(k) - zu(k) )**2 ) )
1616!
1617!--                   xy-distance (horizontal edges)
1618                      IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 )  .OR.        &
1619                           .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.        &
1620                           .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.        &
1621                           .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) )            &
1622                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
1623                                              SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
1624!
1625!--                   xyz distance (vertical and horizontal edges, down)
1626                      IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 )  .OR.      &
1627                           .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.      &
1628                           .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.      &
1629                           .NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) )          &
1630                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
1631                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
1632                                                    +  ( zu(k) - zw(k-1) )**2  ) )
1633!
1634!--                   xyz distance (vertical and horizontal edges, up)
1635                      IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 )  .OR.      &
1636                           .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.      &
1637                           .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.      &
1638                           .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) )          &
1639                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
1640                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
1641                                                    +  ( zw(k) - zu(k) )**2  ) )
1642
1643                   ENDIF
1644!
1645!--                Adjust mixing length by wall-adjustment factor and limit it by l_grid
1646                   l_wall(k,j,i) = MIN( l_wall(k,j,i) * wall_adjustment_factor, l_grid(k) )
1647
1648                ENDDO  !k loop
1649             ENDDO  !j loop
1650          ENDDO  !i loop
1651
1652       ENDIF  !if wall_adjustment
1653
1654    ELSE
1655!
1656!--    Initialize the mixing length in case of a RANS simulation
1657       ALLOCATE( l_black(nzb:nzt+1) )
1658
1659!
1660!--    Calculate mixing length according to Blackadar (1962)
1661       IF ( f /= 0.0_wp )  THEN
1662          l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /            &
1663                  ABS( f ) + 1.0E-10_wp
1664       ELSE
1665          l_max = 30.0_wp
1666       ENDIF
1667
1668       DO  k = nzb, nzt
1669          l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max )
1670       ENDDO
1671
1672       l_black(nzt+1) = l_black(nzt)
1673
1674!
1675!--    Get height level of highest topography within local subdomain
1676       k_max_topo = 0
1677       DO  i = nxlg, nxrg
1678          DO  j = nysg, nyng
1679             DO  k = nzb+1, nzt-1
1680                IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND.  &
1681                     k > k_max_topo )  &
1682                   k_max_topo = k
1683             ENDDO
1684          ENDDO
1685       ENDDO
1686
1687       l_wall(nzb,:,:) = l_black(nzb)
1688       l_wall(nzt+1,:,:) = l_black(nzt+1)
1689!
1690!--    Limit mixing length to either nearest wall or Blackadar mixing length.
1691!--    For that, analyze each grid point (i/j/k) ("analysed grid point") and
1692!--    search within its vicinity for the shortest distance to a wall by cal-
1693!--    culating the distance between the analysed grid point and the "viewed
1694!--    grid point" if it contains a wall (belongs to topography).
1695       DO  k = nzb+1, nzt
1696
1697          radius = l_black(k)  ! radius within walls are searched
1698!
1699!--       Set l_wall to its default maximum value (l_back)
1700          l_wall(k,:,:) = radius
1701
1702!
1703!--       Compute search radius as number of grid points in all directions
1704          rad_i = CEILING( radius / dx )
1705          rad_j = CEILING( radius / dy )
1706
1707          DO  kk = 0, nzt-k
1708             rad_k_t = kk
1709!
1710!--          Limit upward search radius to height of maximum topography
1711             IF ( zu(k+kk)-zu(k) >= radius .OR. k+kk >= k_max_topo )  EXIT
1712          ENDDO
1713
1714          DO  kk = 0, k
1715             rad_k_b = kk
1716             IF ( zu(k)-zu(k-kk) >= radius )  EXIT
1717          ENDDO
1718
1719!
1720!--       Get maximum vertical radius; necessary for defining arrays
1721          rad_k = MAX( rad_k_b, rad_k_t )
1722!
1723!--       When analysed grid point lies above maximum topography, set search
1724!--       radius to 0 if the distance between the analysed grid point and max
1725!--       topography height is larger than the maximum search radius
1726          IF ( zu(k-rad_k_b) > zu(k_max_topo) )  rad_k_b = 0
1727!
1728!--       Search within vicinity only if the vertical search radius is >0
1729          IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 )  THEN
1730
1731             !> @note shape of vicinity is larger in z direction
1732             !>   Shape of vicinity is two grid points larger than actual search
1733             !>   radius in vertical direction. The first and last grid point is
1734             !>   always set to 1 to asure correct detection of topography. See
1735             !>   function "shortest_distance" for details.
1736             !>   2018-03-16, gronemeier
1737             ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
1738             ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
1739
1740             vicinity = 1
1741
1742             DO  i = nxl, nxr
1743                DO  j = nys, nyn
1744!
1745!--                Start search only if (i/j/k) belongs to atmosphere
1746                   IF ( BTEST( wall_flags_0(k,j,i), 0 )  )  THEN
1747!
1748!--                   Reset topography within vicinity
1749                      vicinity(-rad_k:rad_k,:,:) = 0
1750!
1751!--                   Copy area surrounding analysed grid point into vicinity.
1752!--                   First, limit size of data copied to vicinity by the domain
1753!--                   border
1754                      !> @note limit copied area to 1 grid point in hor. dir.
1755                      !>   Ignore walls in horizontal direction which are
1756                      !>   further away than a single grid point. This allows to
1757                      !>   only search within local subdomain without the need
1758                      !>   of global topography information.
1759                      !>   The error made by this assumption are acceptable at
1760                      !>   the moment.
1761                      !>   2018-10-01, gronemeier
1762                      rad_i_l = MIN( 1, rad_i, i )
1763                      rad_i_r = MIN( 1, rad_i, nx-i )
1764
1765                      rad_j_s = MIN( 1, rad_j, j )
1766                      rad_j_n = MIN( 1, rad_j, ny-j )
1767
1768                      CALL copy_into_vicinity( k, j, i,           &
1769                                               -rad_k_b, rad_k_t, &
1770                                               -rad_j_s, rad_j_n, &
1771                                               -rad_i_l, rad_i_r  )
1772                      !> @note in case of cyclic boundaries, those parts of the
1773                      !>   topography which lies beyond the domain borders but
1774                      !>   still within the search radius still needs to be
1775                      !>   copied into 'vicinity'. As the effective search
1776                      !>   radius is limited to 1 at the moment, no further
1777                      !>   copying is needed. Old implementation (prior to
1778                      !>   2018-10-01) had this covered but used a global array.
1779                      !>   2018-10-01, gronemeier
1780
1781!
1782!--                   Search for walls only if there is any within vicinity
1783                      IF ( MAXVAL( vicinity(-rad_k:rad_k,:,:) ) /= 0 )  THEN
1784!
1785!--                      Search within first half (positive x)
1786                         dist_dx = rad_i
1787                         DO  ii = 0, dist_dx
1788!
1789!--                         Search along vertical direction only if below
1790!--                         maximum topography
1791                            IF ( rad_k_t > 0 ) THEN
1792!
1793!--                            Search for walls within octant (+++)
1794                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1795                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1796                                       shortest_distance( vic_yz, .TRUE., ii ) )
1797!
1798!--                            Search for walls within octant (+-+)
1799!--                            Switch order of array so that the analysed grid
1800!--                            point is always located at (0/0) (required by
1801!--                            shortest_distance").
1802                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1803                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1804                                       shortest_distance( vic_yz, .TRUE., ii ) )
1805
1806                            ENDIF
1807!
1808!--                         Search for walls within octant (+--)
1809                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1810                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1811                                      shortest_distance( vic_yz, .FALSE., ii ) )
1812!
1813!--                         Search for walls within octant (++-)
1814                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1815                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1816                                      shortest_distance( vic_yz, .FALSE., ii ) )
1817!
1818!--                         Reduce search along x by already found distance
1819                            dist_dx = CEILING( l_wall(k,j,i) / dx )
1820
1821                         ENDDO
1822!
1823!-                       Search within second half (negative x)
1824                         DO  ii = 0, -dist_dx, -1
1825!
1826!--                         Search along vertical direction only if below
1827!--                         maximum topography
1828                            IF ( rad_k_t > 0 ) THEN
1829!
1830!--                            Search for walls within octant (-++)
1831                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1832                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1833                                      shortest_distance( vic_yz, .TRUE., -ii ) )
1834!
1835!--                            Search for walls within octant (--+)
1836!--                            Switch order of array so that the analysed grid
1837!--                            point is always located at (0/0) (required by
1838!--                            shortest_distance").
1839                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1840                               l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
1841                                      shortest_distance( vic_yz, .TRUE., -ii ) )
1842
1843                            ENDIF
1844!
1845!--                         Search for walls within octant (---)
1846                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1847                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1848                                     shortest_distance( vic_yz, .FALSE., -ii ) )
1849!
1850!--                         Search for walls within octant (-+-)
1851                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1852                            l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
1853                                     shortest_distance( vic_yz, .FALSE., -ii ) )
1854!
1855!--                         Reduce search along x by already found distance
1856                            dist_dx = CEILING( l_wall(k,j,i) / dx )
1857
1858                         ENDDO
1859
1860                      ENDIF  !Check for any walls within vicinity
1861
1862                   ELSE  !Check if (i,j,k) belongs to atmosphere
1863
1864                      l_wall(k,j,i) = l_black(k)
1865
1866                   ENDIF
1867
1868                ENDDO  !j loop
1869             ENDDO  !i loop
1870
1871             DEALLOCATE( vicinity )
1872             DEALLOCATE( vic_yz )
1873
1874          ENDIF  !check vertical size of vicinity
1875
1876       ENDDO  !k loop
1877
1878       !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
1879
1880    ENDIF  !LES or RANS mode
1881
1882!
1883!-- Set lateral boundary conditions for l_wall
1884    CALL exchange_horiz( l_wall, nbgp )
1885
1886    !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
1887    !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
1888
1889    CONTAINS
1890!------------------------------------------------------------------------------!
1891! Description:
1892! ------------
1893!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and
1894!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
1895!> closest to the origin (0/0) of 'array'.
1896!------------------------------------------------------------------------------!
1897    REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
1898
1899       IMPLICIT NONE
1900
1901       LOGICAL, INTENT(IN) :: orientation    !< flag if array represents an array oriented upwards (true) or downwards (false)
1902
1903       INTEGER(iwp), INTENT(IN) :: pos_i     !< x position of the yz-plane 'array'
1904
1905       INTEGER(iwp) :: a                     !< loop index
1906       INTEGER(iwp) :: b                     !< loop index
1907       INTEGER(iwp) :: jj                    !< loop index
1908
1909       INTEGER(KIND=1) :: maximum            !< maximum of array along z dimension
1910
1911       INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
1912
1913       INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i
1914
1915!
1916!--    Get coordinate of first maximum along vertical dimension
1917!--    at each y position of array (similar to function maxloc but more stable).
1918       DO  a = 0, rad_j
1919          loc_k(a) = rad_k+1
1920          maximum = MAXVAL( array(:,a) )
1921          DO  b = 0, rad_k+1
1922             IF ( array(b,a) == maximum )  THEN
1923                loc_k(a) = b
1924                EXIT
1925             ENDIF
1926          ENDDO
1927       ENDDO
1928!
1929!--    Set distance to the default maximum value (=search radius)
1930       shortest_distance = radius
1931!
1932!--    Calculate distance between position (0/0/0) and
1933!--    position (pos_i/jj/loc(jj)) and only save the shortest distance.
1934       IF ( orientation ) THEN  !if array is oriented upwards
1935          DO  jj = 0, rad_j
1936             shortest_distance =                                               &
1937                MIN( shortest_distance,                                        &
1938                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
1939                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
1940                         + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2             &
1941                         )                                                     &
1942                   )
1943          ENDDO
1944       ELSE  !if array is oriented downwards
1945          !> @note MAX within zw required to circumvent error at domain border
1946          !>   At the domain border, if non-cyclic boundary is present, the
1947          !>   index for zw could be -1, which will be errorneous (zw(-1) does
1948          !>   not exist). The MAX function limits the index to be at least 0.
1949          DO  jj = 0, rad_j
1950             shortest_distance =                                               &
1951                MIN( shortest_distance,                                        &
1952                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
1953                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
1954                         + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2    &
1955                         )                                                     &
1956                   )
1957          ENDDO
1958       ENDIF
1959
1960    END FUNCTION
1961
1962!------------------------------------------------------------------------------!
1963! Description:
1964! ------------
1965!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point
1966!> (kp,jp,ip) containing the first bit of wall_flags_0 into the array
1967!> 'vicinity'. Only copy first bit as this indicates the presence of topography.
1968!------------------------------------------------------------------------------!
1969    SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir )
1970
1971       IMPLICIT NONE
1972
1973       INTEGER(iwp), INTENT(IN) :: il !< left loop boundary
1974       INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction
1975       INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary
1976       INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary
1977       INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction
1978       INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary
1979       INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary
1980       INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction
1981       INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary
1982
1983       INTEGER(iwp) :: i   !< loop index
1984       INTEGER(iwp) :: j   !< loop index
1985       INTEGER(iwp) :: k   !< loop index
1986
1987       DO  i = il, ir
1988          DO  j = js, jn
1989             DO  k = kb, kt
1990                vicinity(k,j,i) = MERGE( 0, 1,               &
1991                       BTEST( wall_flags_0(kp+k,jp+j,ip+i), 0 ) )
1992             ENDDO
1993          ENDDO
1994       ENDDO
1995
1996    END SUBROUTINE copy_into_vicinity
1997
1998 END SUBROUTINE tcm_init_mixing_length
1999
2000
2001!------------------------------------------------------------------------------!
2002! Description:
2003! ------------
2004!> Initialize virtual velocities used later in production_e.
2005!------------------------------------------------------------------------------!
2006 SUBROUTINE production_e_init
2007
2008    USE arrays_3d,                                                             &
2009        ONLY:  drho_air_zw, zu
2010
2011    USE control_parameters,                                                    &
2012        ONLY:  constant_flux_layer
2013
2014    USE surface_layer_fluxes_mod,                                              &
2015        ONLY:  phi_m
2016
2017    IMPLICIT NONE
2018
2019    INTEGER(iwp) ::  i      !< grid index x-direction
2020    INTEGER(iwp) ::  j      !< grid index y-direction
2021    INTEGER(iwp) ::  k      !< grid index z-direction
2022    INTEGER(iwp) ::  m      !< running index surface elements
2023
2024    REAL(wp) ::  km_sfc     !< diffusion coefficient, used to compute virtual velocities
2025
2026    IF ( constant_flux_layer )  THEN
2027!
2028!--    Calculate a virtual velocity at the surface in a way that the
2029!--    vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
2030!--    Prandtl law (-w'u'/km). This gradient is used in the TKE shear
2031!--    production term at k=1 (see production_e_ij).
2032!--    The velocity gradient has to be limited in case of too small km
2033!--    (otherwise the timestep may be significantly reduced by large
2034!--    surface winds).
2035!--    not available in case of non-cyclic boundary conditions.
2036!--    Default surfaces, upward-facing
2037       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
2038       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
2039       !$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu)
2040       DO  m = 1, surf_def_h(0)%ns
2041
2042          i = surf_def_h(0)%i(m)
2043          j = surf_def_h(0)%j(m)
2044          k = surf_def_h(0)%k(m)
2045!
2046!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
2047!--       and km are not on the same grid. Actually, a further
2048!--       interpolation of km onto the u/v-grid is necessary. However, the
2049!--       effect of this error is negligible.
2050          km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) /       &
2051                   phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
2052
2053          surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
2054                                     drho_air_zw(k-1)               *          &
2055                                     ( zu(k+1) - zu(k-1)    )       /          &
2056                                     ( km_sfc  + 1.0E-20_wp )
2057          surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
2058                                     drho_air_zw(k-1)               *          &
2059                                     ( zu(k+1) - zu(k-1)    )       /          &
2060                                     ( km_sfc  + 1.0E-20_wp )
2061
2062          IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
2063               ABS( u(k+1,j,i) - u(k-1,j,i)           )                        &
2064             )  surf_def_h(0)%u_0(m) = u(k-1,j,i)
2065
2066          IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) )  >                     &
2067               ABS( v(k+1,j,i) - v(k-1,j,i)           )                        &
2068             )  surf_def_h(0)%v_0(m) = v(k-1,j,i)
2069
2070       ENDDO
2071!
2072!--    Default surfaces, downward-facing surfaces
2073       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
2074       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
2075       !$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km)
2076       DO  m = 1, surf_def_h(1)%ns
2077
2078          i = surf_def_h(1)%i(m)
2079          j = surf_def_h(1)%j(m)
2080          k = surf_def_h(1)%k(m)
2081!
2082!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
2083!--       and km are not on the same grid. Actually, a further
2084!--       interpolation of km onto the u/v-grid is necessary. However, the
2085!--       effect of this error is negligible.
2086          surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *          &
2087                                     drho_air_zw(k-1) *                        &
2088                                     ( zu(k+1)    - zu(k-1)    )  /            &
2089                                     ( km(k,j,i)  + 1.0E-20_wp )
2090          surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *          &
2091                                     drho_air_zw(k-1) *                        &
2092                                     ( zu(k+1)    - zu(k-1)    )  /            &
2093                                     ( km(k,j,i)  + 1.0E-20_wp )
2094
2095          IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                     &
2096               ABS( u(k+1,j,i)           - u(k-1,j,i) )                        &
2097             )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
2098
2099          IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >                     &
2100               ABS( v(k+1,j,i)           - v(k-1,j,i) )                        &
2101             )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
2102
2103       ENDDO
2104!
2105!--    Natural surfaces, upward-facing
2106       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
2107       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
2108       !$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu)
2109       DO  m = 1, surf_lsm_h%ns
2110
2111          i = surf_lsm_h%i(m)
2112          j = surf_lsm_h%j(m)
2113          k = surf_lsm_h%k(m)
2114!
2115!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
2116!--       and km are not on the same grid. Actually, a further
2117!--       interpolation of km onto the u/v-grid is necessary. However, the
2118!--       effect of this error is negligible.
2119          km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /             &
2120                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
2121
2122          surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
2123                                        drho_air_zw(k-1)         *             &
2124                                        ( zu(k+1) - zu(k-1)    ) /             &
2125                                        ( km_sfc  + 1.0E-20_wp )
2126          surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
2127                                        drho_air_zw(k-1)         *             &
2128                                        ( zu(k+1) - zu(k-1)    ) /             &
2129                                        ( km_sfc  + 1.0E-20_wp )
2130
2131          IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
2132               ABS( u(k+1,j,i) - u(k-1,j,i)   )                                &
2133             )  surf_lsm_h%u_0(m) = u(k-1,j,i)
2134
2135          IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) )  >                        &
2136               ABS( v(k+1,j,i) - v(k-1,j,i)   )                                &
2137             )  surf_lsm_h%v_0(m) = v(k-1,j,i)
2138
2139       ENDDO
2140!
2141!--    Urban surfaces, upward-facing
2142       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
2143       !$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
2144       !$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu)
2145       DO  m = 1, surf_usm_h%ns
2146
2147          i = surf_usm_h%i(m)
2148          j = surf_usm_h%j(m)
2149          k = surf_usm_h%k(m)
2150!
2151!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
2152!--       and km are not on the same grid. Actually, a further
2153!--       interpolation of km onto the u/v-grid is necessary. However, the
2154!--       effect of this error is negligible.
2155          km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /             &
2156                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
2157
2158          surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)    *             &
2159                                        drho_air_zw(k-1)         *             &
2160                                        ( zu(k+1) - zu(k-1)    ) /             &
2161                                        ( km_sfc  + 1.0E-20_wp )
2162          surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)    *             &
2163                                        drho_air_zw(k-1)         *             &
2164                                        ( zu(k+1) - zu(k-1)    ) /             &
2165                                        ( km_sfc  + 1.0E-20_wp )
2166
2167          IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
2168               ABS( u(k+1,j,i) - u(k-1,j,i)   )                                &
2169             )  surf_usm_h%u_0(m) = u(k-1,j,i)
2170
2171          IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) )  >                        &
2172               ABS( v(k+1,j,i) - v(k-1,j,i)   )                                &
2173             )  surf_usm_h%v_0(m) = v(k-1,j,i)
2174
2175       ENDDO
2176
2177    ENDIF
2178
2179 END SUBROUTINE production_e_init
2180
2181
2182!--------------------------------------------------------------------------------------------------!
2183! Description:
2184! ------------
2185!> Execute module-specific actions for all grid points
2186!--------------------------------------------------------------------------------------------------!
2187 SUBROUTINE tcm_actions( location )
2188
2189
2190    CHARACTER (LEN=*) ::  location !<
2191
2192!    INTEGER(iwp) ::  i !<
2193!    INTEGER(iwp) ::  j !<
2194!    INTEGER(iwp) ::  k !<
2195
2196!
2197!-- Here the module-specific actions follow
2198!-- No calls for single grid points are allowed at locations before and
2199!-- after the timestep, since these calls are not within an i,j-loop
2200    SELECT CASE ( location )
2201
2202       CASE ( 'before_timestep' )
2203
2204
2205       CASE ( 'before_prognostic_equations' )
2206
2207          IF ( .NOT. constant_diffusion )  CALL production_e_init
2208
2209
2210       CASE ( 'after_integration' )
2211
2212
2213       CASE ( 'after_timestep' )
2214
2215
2216       CASE ( 'u-tendency' )
2217
2218
2219       CASE ( 'v-tendency' )
2220
2221
2222       CASE ( 'w-tendency' )
2223
2224
2225       CASE ( 'pt-tendency' )
2226
2227
2228       CASE ( 'sa-tendency' )
2229
2230
2231       CASE ( 'e-tendency' )
2232
2233
2234       CASE ( 'q-tendency' )
2235
2236
2237       CASE ( 's-tendency' )
2238
2239
2240       CASE DEFAULT
2241          CONTINUE
2242
2243    END SELECT
2244
2245 END SUBROUTINE tcm_actions
2246
2247
2248!--------------------------------------------------------------------------------------------------!
2249! Description:
2250! ------------
2251!> Execute module-specific actions for grid point i,j
2252!--------------------------------------------------------------------------------------------------!
2253 SUBROUTINE tcm_actions_ij( i, j, location )
2254
2255
2256    CHARACTER (LEN=*) ::  location
2257
2258    INTEGER(iwp) ::  i
2259    INTEGER(iwp) ::  j
2260
2261!
2262!-- Here the module-specific actions follow
2263    SELECT CASE ( location )
2264
2265       CASE ( 'u-tendency' )
2266
2267!--       Next line is to avoid compiler warning about unused variables. Please remove.
2268          IF ( i +  j < 0 )  CONTINUE
2269
2270       CASE ( 'v-tendency' )
2271
2272
2273       CASE ( 'w-tendency' )
2274
2275
2276       CASE ( 'pt-tendency' )
2277
2278
2279       CASE ( 'sa-tendency' )
2280
2281
2282       CASE ( 'e-tendency' )
2283
2284
2285       CASE ( 'q-tendency' )
2286
2287
2288       CASE ( 's-tendency' )
2289
2290
2291       CASE DEFAULT
2292          CONTINUE
2293
2294    END SELECT
2295
2296 END SUBROUTINE tcm_actions_ij
2297
2298
2299!------------------------------------------------------------------------------!
2300! Description:
2301! ------------
2302!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2303!> Vector-optimized version
2304!------------------------------------------------------------------------------!
2305 SUBROUTINE tcm_prognostic_equations
2306
2307    USE control_parameters,                                                    &
2308        ONLY:  scalar_advec, tsc
2309
2310    IMPLICIT NONE
2311
2312    INTEGER(iwp) ::  i       !< loop index
2313    INTEGER(iwp) ::  j       !< loop index
2314    INTEGER(iwp) ::  k       !< loop index
2315
2316    REAL(wp)     ::  sbt     !< wheighting factor for sub-time step
2317
2318!
2319!-- If required, compute prognostic equation for turbulent kinetic
2320!-- energy (TKE)
2321    IF ( .NOT. constant_diffusion )  THEN
2322
2323       CALL cpu_log( log_point_s(67), 'tke-equation', 'start' )
2324
2325       sbt = tsc(2)
2326       IF ( .NOT. use_upstream_for_tke )  THEN
2327          IF ( scalar_advec == 'bc-scheme' )  THEN
2328
2329             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2330!
2331!--             Bott-Chlond scheme always uses Euler time step. Thus:
2332                sbt = 1.0_wp
2333             ENDIF
2334             tend = 0.0_wp
2335             CALL advec_s_bc( e, 'e' )
2336
2337          ENDIF
2338       ENDIF
2339
2340!
2341!--    TKE-tendency terms with no communication
2342       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
2343          IF ( use_upstream_for_tke )  THEN
2344             tend = 0.0_wp
2345             CALL advec_s_up( e )
2346          ELSE
2347             !$ACC KERNELS PRESENT(tend)
2348             tend = 0.0_wp
2349             !$ACC END KERNELS
2350             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2351                IF ( ws_scheme_sca )  THEN
2352                   CALL advec_s_ws( advc_flags_s, e, 'e',                      &
2353                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
2354                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
2355                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
2356                                    bc_dirichlet_s  .OR.  bc_radiation_s )
2357                ELSE
2358                   CALL advec_s_pw( e )
2359                ENDIF
2360             ELSE
2361                CALL advec_s_up( e )
2362             ENDIF
2363          ENDIF
2364       ENDIF
2365
2366       CALL production_e( .FALSE. )
2367
2368       IF ( .NOT. humidity )  THEN
2369          IF ( ocean_mode )  THEN
2370             CALL diffusion_e( prho, prho_reference )
2371          ELSE
2372             CALL diffusion_e( pt, pt_reference )
2373          ENDIF
2374       ELSE
2375          CALL diffusion_e( vpt, pt_reference )
2376       ENDIF
2377
2378!
2379!--    Additional sink term for flows through plant canopies
2380       IF ( plant_canopy )  CALL pcm_tendency( 6 )
2381
2382!       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
2383
2384!
2385!--    Prognostic equation for TKE.
2386!--    Eliminate negative TKE values, which can occur due to numerical
2387!--    reasons in the course of the integration. In such cases the old TKE
2388!--    value is reduced by 90%.
2389       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2390       !$ACC PRESENT(e, tend, te_m, wall_flags_0) &
2391       !$ACC PRESENT(tsc(3:3)) &
2392       !$ACC PRESENT(e_p)
2393       DO  i = nxl, nxr
2394          DO  j = nys, nyn
2395             DO  k = nzb+1, nzt
2396                e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
2397                                                 tsc(3) * te_m(k,j,i) )        &
2398                                        )                                      &
2399                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2400                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2401                                          )
2402                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
2403             ENDDO
2404          ENDDO
2405       ENDDO
2406
2407!
2408!--    Calculate tendencies for the next Runge-Kutta step
2409       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2410          IF ( intermediate_timestep_count == 1 )  THEN
2411             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2412             !$ACC PRESENT(tend, te_m)
2413             DO  i = nxl, nxr
2414                DO  j = nys, nyn
2415                   DO  k = nzb+1, nzt
2416                      te_m(k,j,i) = tend(k,j,i)
2417                   ENDDO
2418                ENDDO
2419             ENDDO
2420          ELSEIF ( intermediate_timestep_count < &
2421                   intermediate_timestep_count_max )  THEN
2422             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2423             !$ACC PRESENT(tend, te_m)
2424             DO  i = nxl, nxr
2425                DO  j = nys, nyn
2426                   DO  k = nzb+1, nzt
2427                      te_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
2428                                     + 5.3125_wp * te_m(k,j,i)
2429                   ENDDO
2430                ENDDO
2431             ENDDO
2432          ENDIF
2433       ENDIF
2434
2435       CALL cpu_log( log_point_s(67), 'tke-equation', 'stop' )
2436
2437    ENDIF   ! TKE equation
2438
2439!
2440!-- If required, compute prognostic equation for TKE dissipation rate
2441    IF ( rans_tke_e )  THEN
2442
2443       CALL cpu_log( log_point_s(64), 'diss-equation', 'start' )
2444
2445       sbt = tsc(2)
2446       IF ( .NOT. use_upstream_for_tke )  THEN
2447          IF ( scalar_advec == 'bc-scheme' )  THEN
2448
2449             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
2450!
2451!--             Bott-Chlond scheme always uses Euler time step. Thus:
2452                sbt = 1.0_wp
2453             ENDIF
2454             tend = 0.0_wp
2455             CALL advec_s_bc( diss, 'diss' )
2456
2457          ENDIF
2458       ENDIF
2459
2460!
2461!--    dissipation-tendency terms with no communication
2462       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
2463          IF ( use_upstream_for_tke )  THEN
2464             tend = 0.0_wp
2465             CALL advec_s_up( diss )
2466          ELSE
2467             tend = 0.0_wp
2468             IF ( timestep_scheme(1:5) == 'runge' )  THEN
2469                IF ( ws_scheme_sca )  THEN
2470                   CALL advec_s_ws( advc_flags_s, diss, 'diss',                &
2471                                    bc_dirichlet_l  .OR.  bc_radiation_l,      &
2472                                    bc_dirichlet_n  .OR.  bc_radiation_n,      &
2473                                    bc_dirichlet_r  .OR.  bc_radiation_r,      &
2474                                    bc_dirichlet_s  .OR.  bc_radiation_s )
2475                ELSE
2476                   CALL advec_s_pw( diss )
2477                ENDIF
2478             ELSE
2479                CALL advec_s_up( diss )
2480             ENDIF
2481          ENDIF
2482       ENDIF
2483!
2484!--    Production of TKE dissipation rate
2485       CALL production_e( .TRUE. )
2486!
2487!--    Diffusion term of TKE dissipation rate
2488       CALL diffusion_diss
2489!
2490!--    Additional sink term for flows through plant canopies
2491!        IF ( plant_canopy )  CALL pcm_tendency( ? )         !> @todo not yet implemented
2492
2493!       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
2494
2495!
2496!--    Prognostic equation for TKE dissipation.
2497!--    Eliminate negative dissipation values, which can occur due to numerical
2498!--    reasons in the course of the integration. In such cases the old
2499!--    dissipation value is reduced by 90%.
2500       DO  i = nxl, nxr
2501          DO  j = nys, nyn
2502             DO  k = nzb+1, nzt
2503                diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +  &
2504                                                 tsc(3) * tdiss_m(k,j,i) )     &
2505                                        )                                      &
2506                                   * MERGE( 1.0_wp, 0.0_wp,                    &
2507                                             BTEST( wall_flags_0(k,j,i), 0 )   &
2508                                          )
2509                IF ( diss_p(k,j,i) < 0.0_wp )                                  &
2510                   diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
2511             ENDDO
2512          ENDDO
2513       ENDDO
2514
2515!
2516!--    Calculate tendencies for the next Runge-Kutta step
2517       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2518          IF ( intermediate_timestep_count == 1 )  THEN
2519             DO  i = nxl, nxr
2520                DO  j = nys, nyn
2521                   DO  k = nzb+1, nzt
2522                      tdiss_m(k,j,i) = tend(k,j,i)
2523                   ENDDO
2524                ENDDO
2525             ENDDO
2526          ELSEIF ( intermediate_timestep_count < &
2527                   intermediate_timestep_count_max )  THEN
2528             DO  i = nxl, nxr
2529                DO  j = nys, nyn
2530                   DO  k = nzb+1, nzt
2531                      tdiss_m(k,j,i) =   -9.5625_wp * tend(k,j,i)              &
2532                                        + 5.3125_wp * tdiss_m(k,j,i)
2533                   ENDDO
2534                ENDDO
2535             ENDDO
2536          ENDIF
2537       ENDIF
2538
2539       CALL cpu_log( log_point_s(64), 'diss-equation', 'stop' )
2540
2541    ENDIF
2542
2543 END SUBROUTINE tcm_prognostic_equations
2544
2545
2546!------------------------------------------------------------------------------!
2547! Description:
2548! ------------
2549!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2550!> Cache-optimized version
2551!------------------------------------------------------------------------------!
2552 SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn )
2553
2554    USE arrays_3d,                                                             &
2555        ONLY:  diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss,      &
2556               flux_l_e, flux_s_diss, flux_s_e
2557
2558    USE control_parameters,                                                    &
2559        ONLY:  tsc
2560
2561    IMPLICIT NONE
2562
2563    INTEGER(iwp) ::  i       !< loop index x direction
2564    INTEGER(iwp) ::  i_omp   !< first loop index of i-loop in prognostic_equations
2565    INTEGER(iwp) ::  j       !< loop index y direction
2566    INTEGER(iwp) ::  k       !< loop index z direction
2567    INTEGER(iwp) ::  tn      !< task number of openmp task
2568
2569!
2570!-- If required, compute prognostic equation for turbulent kinetic
2571!-- energy (TKE)
2572    IF ( .NOT. constant_diffusion )  THEN
2573
2574!
2575!--    Tendency-terms for TKE
2576       tend(:,j,i) = 0.0_wp
2577       IF ( timestep_scheme(1:5) == 'runge'  &
2578           .AND.  .NOT. use_upstream_for_tke )  THEN
2579           IF ( ws_scheme_sca )  THEN
2580               CALL advec_s_ws( advc_flags_s,                                  &
2581                                i, j, e, 'e', flux_s_e, diss_s_e,              &
2582                                flux_l_e, diss_l_e , i_omp, tn,                &
2583                                bc_dirichlet_l  .OR.  bc_radiation_l,          &
2584                                bc_dirichlet_n  .OR.  bc_radiation_n,          &
2585                                bc_dirichlet_r  .OR.  bc_radiation_r,          &
2586                                bc_dirichlet_s  .OR.  bc_radiation_s )
2587           ELSE
2588               CALL advec_s_pw( i, j, e )
2589           ENDIF
2590       ELSE
2591          CALL advec_s_up( i, j, e )
2592       ENDIF
2593
2594       CALL production_e_ij( i, j, .FALSE. )
2595
2596       IF ( .NOT. humidity )  THEN
2597          IF ( ocean_mode )  THEN
2598             CALL diffusion_e_ij( i, j, prho, prho_reference )
2599          ELSE
2600             CALL diffusion_e_ij( i, j, pt, pt_reference )
2601          ENDIF
2602       ELSE
2603          CALL diffusion_e_ij( i, j, vpt, pt_reference )
2604       ENDIF
2605
2606!
2607!--    Additional sink term for flows through plant canopies
2608       IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
2609
2610!       CALL user_actions( i, j, 'e-tendency' ) ToDo: find general solution for circular dependency between modules
2611
2612!
2613!--    Prognostic equation for TKE.
2614!--    Eliminate negative TKE values, which can occur due to numerical
2615!--    reasons in the course of the integration. In such cases the old
2616!--    TKE value is reduced by 90%.
2617       DO  k = nzb+1, nzt
2618          e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +           &
2619                                              tsc(3) * te_m(k,j,i) )           &
2620                                  )                                            &
2621                                 * MERGE( 1.0_wp, 0.0_wp,                      &
2622                                          BTEST( wall_flags_0(k,j,i), 0 )      &
2623                                        )
2624          IF ( e_p(k,j,i) <= 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
2625       ENDDO
2626
2627!
2628!--    Calculate tendencies for the next Runge-Kutta step
2629       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2630          IF ( intermediate_timestep_count == 1 )  THEN
2631             DO  k = nzb+1, nzt
2632                te_m(k,j,i) = tend(k,j,i)
2633             ENDDO
2634          ELSEIF ( intermediate_timestep_count < &
2635                   intermediate_timestep_count_max )  THEN
2636             DO  k = nzb+1, nzt
2637                te_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                     &
2638                                 5.3125_wp * te_m(k,j,i)
2639             ENDDO
2640          ENDIF
2641       ENDIF
2642
2643    ENDIF   ! TKE equation
2644
2645!
2646!-- If required, compute prognostic equation for TKE dissipation rate
2647    IF ( rans_tke_e )  THEN
2648!
2649!--    Tendency-terms for dissipation
2650       tend(:,j,i) = 0.0_wp
2651       IF ( timestep_scheme(1:5) == 'runge'  &
2652           .AND.  .NOT. use_upstream_for_tke )  THEN
2653           IF ( ws_scheme_sca )  THEN
2654               CALL advec_s_ws( advc_flags_s,                                  &
2655                                i, j, diss, 'diss', flux_s_diss, diss_s_diss,  &
2656                                flux_l_diss, diss_l_diss, i_omp, tn,           &
2657                                bc_dirichlet_l  .OR.  bc_radiation_l,          &
2658                                bc_dirichlet_n  .OR.  bc_radiation_n,          &
2659                                bc_dirichlet_r  .OR.  bc_radiation_r,          &
2660                                bc_dirichlet_s  .OR.  bc_radiation_s )
2661           ELSE
2662               CALL advec_s_pw( i, j, diss )
2663           ENDIF
2664       ELSE
2665          CALL advec_s_up( i, j, diss )
2666       ENDIF
2667!
2668!--    Production of TKE dissipation rate
2669       CALL production_e_ij( i, j, .TRUE. )
2670!
2671!--    Diffusion term of TKE dissipation rate
2672       CALL diffusion_diss_ij( i, j )
2673!
2674!--    Additional sink term for flows through plant canopies
2675!        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )     !> @todo not yet implemented
2676
2677!       CALL user_actions( i, j, 'diss-tendency' ) ToDo: find general solution for circular dependency between modules
2678
2679!
2680!--    Prognostic equation for TKE dissipation
2681!--    Eliminate negative dissipation values, which can occur due to
2682!--    numerical reasons in the course of the integration. In such cases
2683!--    the old dissipation value is reduced by 90%.
2684       DO  k = nzb+1, nzt
2685          diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +     &
2686                                                    tsc(3) * tdiss_m(k,j,i) )  &
2687                                        )                                      &
2688                                        * MERGE( 1.0_wp, 0.0_wp,               &
2689                                                BTEST( wall_flags_0(k,j,i), 0 )&
2690                                               )
2691       ENDDO
2692
2693!
2694!--    Calculate tendencies for the next Runge-Kutta step
2695       IF ( timestep_scheme(1:5) == 'runge' )  THEN
2696          IF ( intermediate_timestep_count == 1 )  THEN
2697             DO  k = nzb+1, nzt
2698                tdiss_m(k,j,i) = tend(k,j,i)
2699             ENDDO
2700          ELSEIF ( intermediate_timestep_count < &
2701                   intermediate_timestep_count_max )  THEN
2702             DO  k = nzb+1, nzt
2703                tdiss_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                  &
2704                                    5.3125_wp * tdiss_m(k,j,i)
2705             ENDDO
2706          ENDIF
2707       ENDIF
2708
2709    ENDIF   ! dissipation equation
2710
2711 END SUBROUTINE tcm_prognostic_equations_ij
2712
2713
2714!------------------------------------------------------------------------------!
2715! Description:
2716! ------------
2717!> Production terms (shear + buoyancy) of the TKE.
2718!> Vector-optimized version
2719!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
2720!>          not considered well!
2721!------------------------------------------------------------------------------!
2722 SUBROUTINE production_e( diss_production )
2723
2724    USE arrays_3d,                                                             &
2725        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
2726
2727    USE control_parameters,                                                    &
2728        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
2729               rho_reference, use_single_reference_value, use_surface_fluxes,  &
2730               use_top_fluxes
2731
2732    USE grid_variables,                                                        &
2733        ONLY:  ddx, dx, ddy, dy
2734
2735    USE bulk_cloud_model_mod,                                                  &
2736        ONLY:  bulk_cloud_model
2737
2738    IMPLICIT NONE
2739
2740    LOGICAL :: diss_production
2741
2742    INTEGER(iwp) ::  i       !< running index x-direction
2743    INTEGER(iwp) ::  j       !< running index y-direction
2744    INTEGER(iwp) ::  k       !< running index z-direction
2745    INTEGER(iwp) ::  l       !< running index for different surface type orientation
2746    INTEGER(iwp) ::  m       !< running index surface elements
2747    INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
2748    INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
2749    INTEGER(iwp) ::  flag_nr !< number of masking flag
2750
2751    REAL(wp)     ::  def         !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
2752    REAL(wp)     ::  flag        !< flag to mask topography
2753    REAL(wp)     ::  k1          !< temporary factor
2754    REAL(wp)     ::  k2          !< temporary factor
2755    REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
2756    REAL(wp)     ::  theta       !< virtual potential temperature
2757    REAL(wp)     ::  temp        !< theta * Exner-function
2758    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
2759    REAL(wp)     ::  usvs        !< momentum flux u"v"
2760    REAL(wp)     ::  vsus        !< momentum flux v"u"
2761    REAL(wp)     ::  wsus        !< momentum flux w"u"
2762    REAL(wp)     ::  wsvs        !< momentum flux w"v"
2763
2764    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
2765    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
2766    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
2767    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
2768    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
2769    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
2770    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
2771    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
2772    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
2773    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
2774
2775
2776
2777!
2778!-- Calculate TKE production by shear. Calculate gradients at all grid
2779!-- points first, gradients at surface-bounded grid points will be
2780!-- overwritten further below.
2781    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) &
2782    !$ACC PRIVATE(surf_s, surf_e) &
2783    !$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) &
2784    !$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_0) &
2785    !$ACC PRESENT(tend) &
2786    !$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) &
2787    !$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
2788    !$ACC PRESENT(surf_usm_h, surf_usm_v(0:3))
2789    DO  i = nxl, nxr
2790       DO  j = nys, nyn
2791          !$ACC LOOP PRIVATE(k)
2792          DO  k = nzb+1, nzt
2793
2794             dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
2795             dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
2796                                   u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
2797             dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
2798                                   u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
2799
2800             dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
2801                                   v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
2802             dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
2803             dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
2804                                     v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
2805
2806             dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
2807                                   w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
2808             dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
2809                                   w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
2810             dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
2811
2812          ENDDO
2813
2814
2815          flag_nr = 29
2816
2817
2818          IF ( constant_flux_layer )  THEN
2819!
2820
2821             flag_nr = 0
2822
2823!--          Position beneath wall
2824!--          (2) - Will allways be executed.
2825!--          'bottom and wall: use u_0,v_0 and wall functions'
2826!
2827!--          Compute gradients at north- and south-facing surfaces.
2828!--          First, for default surfaces, then for urban surfaces.
2829!--          Note, so far no natural vertical surfaces implemented
2830             DO  l = 0, 1
2831                surf_s = surf_def_v(l)%start_index(j,i)
2832                surf_e = surf_def_v(l)%end_index(j,i)
2833                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2834                DO  m = surf_s, surf_e
2835                   k           = surf_def_v(l)%k(m)
2836                   usvs        = surf_def_v(l)%mom_flux_tke(0,m)
2837                   wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
2838
2839                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
2840                                   * 0.5_wp * dy
2841!
2842!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2843                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2844                                     BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2845                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2846                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2847                ENDDO
2848!
2849!--             Natural surfaces
2850                surf_s = surf_lsm_v(l)%start_index(j,i)
2851                surf_e = surf_lsm_v(l)%end_index(j,i)
2852                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2853                DO  m = surf_s, surf_e
2854                   k           = surf_lsm_v(l)%k(m)
2855                   usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
2856                   wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
2857
2858                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
2859                                   * 0.5_wp * dy
2860!
2861!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2862                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2863                                     BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2864                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2865                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2866                ENDDO
2867!
2868!--             Urban surfaces
2869                surf_s = surf_usm_v(l)%start_index(j,i)
2870                surf_e = surf_usm_v(l)%end_index(j,i)
2871                !$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2872                DO  m = surf_s, surf_e
2873                   k           = surf_usm_v(l)%k(m)
2874                   usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
2875                   wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
2876
2877                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
2878                                   * 0.5_wp * dy
2879!
2880!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2881                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2882                                     BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2883                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2884                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2885                ENDDO
2886             ENDDO
2887!
2888!--          Compute gradients at east- and west-facing walls
2889             DO  l = 2, 3
2890                surf_s = surf_def_v(l)%start_index(j,i)
2891                surf_e = surf_def_v(l)%end_index(j,i)
2892                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2893                DO  m = surf_s, surf_e
2894                   k     = surf_def_v(l)%k(m)
2895                   vsus  = surf_def_v(l)%mom_flux_tke(0,m)
2896                   wsus  = surf_def_v(l)%mom_flux_tke(1,m)
2897
2898                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
2899                                      * 0.5_wp * dx
2900!
2901!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2902                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2903                                     BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2904                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2905                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2906                ENDDO
2907!
2908!--             Natural surfaces
2909                surf_s = surf_lsm_v(l)%start_index(j,i)
2910                surf_e = surf_lsm_v(l)%end_index(j,i)
2911                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2912                DO  m = surf_s, surf_e
2913                   k     = surf_lsm_v(l)%k(m)
2914                   vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
2915                   wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
2916
2917                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
2918                                      * 0.5_wp * dx
2919!
2920!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2921                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2922                                     BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2923                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2924                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2925                ENDDO
2926!
2927!--             Urban surfaces
2928                surf_s = surf_usm_v(l)%start_index(j,i)
2929                surf_e = surf_usm_v(l)%end_index(j,i)
2930                !$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2931                DO  m = surf_s, surf_e
2932                   k     = surf_usm_v(l)%k(m)
2933                   vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
2934                   wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
2935
2936                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
2937                                      * 0.5_wp * dx
2938!
2939!--                -1.0 for right-facing wall, 1.0 for left-facing wall
2940                   sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
2941                                     BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2942                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2943                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2944                ENDDO
2945             ENDDO
2946!
2947!--          Compute gradients at upward-facing surfaces
2948             surf_s = surf_def_h(0)%start_index(j,i)
2949             surf_e = surf_def_h(0)%end_index(j,i)
2950             !$ACC LOOP PRIVATE(m, k)
2951             DO  m = surf_s, surf_e
2952                k = surf_def_h(0)%k(m)
2953!
2954!--             Please note, actually, an interpolation of u_0 and v_0
2955!--             onto the grid center would be required. However, this
2956!--             would require several data transfers between 2D-grid and
2957!--             wall type. The effect of this missing interpolation is
2958!--             negligible. (See also production_e_init).
2959                dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
2960                dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
2961
2962             ENDDO
2963!
2964!--          Natural surfaces
2965             surf_s = surf_lsm_h%start_index(j,i)
2966             surf_e = surf_lsm_h%end_index(j,i)
2967             !$ACC LOOP PRIVATE(m, k)
2968             DO  m = surf_s, surf_e
2969                k = surf_lsm_h%k(m)
2970
2971                dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
2972                dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
2973
2974             ENDDO
2975!
2976!--          Urban surfaces
2977             surf_s = surf_usm_h%start_index(j,i)
2978             surf_e = surf_usm_h%end_index(j,i)
2979             !$ACC LOOP PRIVATE(m, k)
2980             DO  m = surf_s, surf_e
2981                k = surf_usm_h%k(m)
2982
2983                dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
2984                dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
2985
2986             ENDDO
2987!
2988!--          Compute gradients at downward-facing walls, only for
2989!--          non-natural default surfaces
2990             surf_s = surf_def_h(1)%start_index(j,i)
2991             surf_e = surf_def_h(1)%end_index(j,i)
2992             !$ACC LOOP PRIVATE(m, k)
2993             DO  m = surf_s, surf_e
2994                k = surf_def_h(1)%k(m)
2995
2996                dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
2997                dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
2998
2999             ENDDO
3000
3001
3002          ENDIF
3003
3004
3005          !$ACC LOOP PRIVATE(k, def, flag)
3006          DO  k = nzb+1, nzt
3007
3008             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
3009                              dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
3010                              dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
3011                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
3012                              dwdy(k)*dvdz(k) )
3013
3014             IF ( def < 0.0_wp )  def = 0.0_wp
3015
3016             flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
3017
3018             IF ( .NOT. diss_production )  THEN
3019
3020!--             Compute tendency for TKE-production from shear
3021                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
3022
3023             ELSE
3024
3025!--             RANS mode: Compute tendency for dissipation-rate-production from shear
3026                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
3027                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
3028
3029             ENDIF
3030
3031          ENDDO
3032
3033
3034       ENDDO
3035    ENDDO
3036
3037!
3038!-- If required, calculate TKE production by buoyancy
3039    IF ( .NOT. neutral )  THEN
3040
3041       IF ( .NOT. humidity )  THEN
3042
3043          IF ( ocean_mode )  THEN
3044!
3045!--          So far in the ocean no special treatment of density flux
3046!--          in the bottom and top surface layer
3047             DO  i = nxl, nxr
3048                DO  j = nys, nyn
3049
3050                   DO  k = nzb+1, nzt
3051                      tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
3052                   ENDDO
3053!
3054!--                Treatment of near-surface grid points, at up- and down-
3055!--                ward facing surfaces
3056                   IF ( use_surface_fluxes )  THEN
3057                      DO  l = 0, 1
3058                         surf_s = surf_def_h(l)%start_index(j,i)
3059                         surf_e = surf_def_h(l)%end_index(j,i)
3060                         DO  m = surf_s, surf_e
3061                            k = surf_def_h(l)%k(m)
3062                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3063                         ENDDO
3064                      ENDDO
3065                   ENDIF
3066
3067                   IF ( use_top_fluxes )  THEN
3068                      surf_s = surf_def_h(2)%start_index(j,i)
3069                      surf_e = surf_def_h(2)%end_index(j,i)
3070                      DO  m = surf_s, surf_e
3071                         k = surf_def_h(2)%k(m)
3072                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3073                      ENDDO
3074                   ENDIF
3075
3076                   IF ( .NOT. diss_production )  THEN
3077
3078!--                   Compute tendency for TKE-production from shear
3079                      DO  k = nzb+1, nzt
3080                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3081                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3082                                       MERGE( rho_reference, prho(k,j,i),       &
3083                                              use_single_reference_value ) )
3084                      ENDDO
3085
3086                   ELSE
3087
3088!--                   RANS mode: Compute tendency for dissipation-rate-production from shear
3089                      DO  k = nzb+1, nzt
3090                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3091                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3092                                       MERGE( rho_reference, prho(k,j,i),       &
3093                                              use_single_reference_value ) ) *  &
3094                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3095                                       c_3
3096                      ENDDO
3097
3098                   ENDIF
3099
3100                ENDDO
3101             ENDDO
3102
3103          ELSE ! or IF ( .NOT. ocean_mode )  THEN
3104
3105             !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
3106             !$ACC PRIVATE(surf_s, surf_e) &
3107             !$ACC PRIVATE(tmp_flux(nzb+1:nzt)) &
3108             !$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_0) &
3109             !$ACC PRESENT(tend) &
3110             !$ACC PRESENT(surf_def_h(0:2)) &
3111             !$ACC PRESENT(surf_lsm_h) &
3112             !$ACC PRESENT(surf_usm_h)
3113             DO  i = nxl, nxr
3114                DO  j = nys, nyn
3115
3116                   !$ACC LOOP PRIVATE(k)
3117                   DO  k = nzb+1, nzt
3118                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
3119                   ENDDO
3120
3121                   IF ( use_surface_fluxes )  THEN
3122!
3123!--                   Default surfaces, up- and downward-facing
3124                      DO  l = 0, 1
3125                         surf_s = surf_def_h(l)%start_index(j,i)
3126                         surf_e = surf_def_h(l)%end_index(j,i)
3127                         !$ACC LOOP PRIVATE(m, k)
3128                         DO  m = surf_s, surf_e
3129                            k = surf_def_h(l)%k(m)
3130                            tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3131                         ENDDO
3132                      ENDDO
3133!
3134!--                   Natural surfaces
3135                      surf_s = surf_lsm_h%start_index(j,i)
3136                      surf_e = surf_lsm_h%end_index(j,i)
3137                      !$ACC LOOP PRIVATE(m, k)
3138                      DO  m = surf_s, surf_e
3139                         k = surf_lsm_h%k(m)
3140                         tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
3141                      ENDDO
3142!
3143!--                   Urban surfaces
3144                      surf_s = surf_usm_h%start_index(j,i)
3145                      surf_e = surf_usm_h%end_index(j,i)
3146                      !$ACC LOOP PRIVATE(m, k)
3147                      DO  m = surf_s, surf_e
3148                         k = surf_usm_h%k(m)
3149                         tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
3150                      ENDDO
3151                   ENDIF
3152
3153                   IF ( use_top_fluxes )  THEN
3154                      surf_s = surf_def_h(2)%start_index(j,i)
3155                      surf_e = surf_def_h(2)%end_index(j,i)
3156                      !$ACC LOOP PRIVATE(m, k)
3157                      DO  m = surf_s, surf_e
3158                         k = surf_def_h(2)%k(m)
3159                         tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3160                      ENDDO
3161                   ENDIF
3162
3163                   IF ( .NOT. diss_production )  THEN
3164
3165!--                   Compute tendency for TKE-production from shear
3166                     !$ACC LOOP PRIVATE(k, flag)
3167                      DO  k = nzb+1, nzt
3168                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3169                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3170                                       MERGE( pt_reference, pt(k,j,i),          &
3171                                              use_single_reference_value ) )
3172                      ENDDO
3173
3174                   ELSE
3175
3176!--                   RANS mode: Compute tendency for dissipation-rate-production from shear
3177                      DO  k = nzb+1, nzt
3178                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3179                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3180                                       MERGE( pt_reference, pt(k,j,i),          &
3181                                              use_single_reference_value ) ) *  &
3182                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3183                                       c_3
3184                      ENDDO
3185
3186                   ENDIF
3187
3188                ENDDO
3189             ENDDO
3190
3191          ENDIF ! from IF ( .NOT. ocean_mode )
3192
3193       ELSE ! or IF ( humidity )  THEN
3194
3195          DO  i = nxl, nxr
3196             DO  j = nys, nyn
3197
3198                DO  k = nzb+1, nzt
3199
3200                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3201                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3202                      k2 = 0.61_wp * pt(k,j,i)
3203                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3204                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3205                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3206                                      ) * dd2zu(k)
3207                   ELSE IF ( bulk_cloud_model )  THEN
3208                      IF ( ql(k,j,i) == 0.0_wp )  THEN
3209                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3210                         k2 = 0.61_wp * pt(k,j,i)
3211                      ELSE
3212                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3213                         temp  = theta * exner(k)
3214                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
3215                                       ( q(k,j,i) - ql(k,j,i) ) *              &
3216                              ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3217                              ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3218                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3219                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3220                      ENDIF
3221                      tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3222                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3223                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3224                                      ) * dd2zu(k)
3225                   ELSE IF ( cloud_droplets )  THEN
3226                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3227                      k2 = 0.61_wp * pt(k,j,i)
3228                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3229                                      ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3230                                        k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
3231                                        pt(k,j,i) * ( ql(k+1,j,i) -            &
3232                                        ql(k-1,j,i) ) ) * dd2zu(k)
3233                   ENDIF
3234
3235                ENDDO
3236
3237                IF ( use_surface_fluxes )  THEN
3238!
3239!--                Treat horizontal default surfaces
3240                   DO  l = 0, 1
3241                      surf_s = surf_def_h(l)%start_index(j,i)
3242                      surf_e = surf_def_h(l)%end_index(j,i)
3243                      DO  m = surf_s, surf_e
3244                         k = surf_def_h(l)%k(m)
3245
3246                         IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3247                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3248                            k2 = 0.61_wp * pt(k,j,i)
3249                         ELSE IF ( bulk_cloud_model )  THEN
3250                            IF ( ql(k,j,i) == 0.0_wp )  THEN
3251                               k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3252                               k2 = 0.61_wp * pt(k,j,i)
3253                            ELSE
3254                               theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3255                               temp  = theta * exner(k)
3256                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
3257                                          ( q(k,j,i) - ql(k,j,i) ) *           &
3258                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3259                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3260                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3261                               k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3262                            ENDIF
3263                         ELSE IF ( cloud_droplets )  THEN
3264                            k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3265                            k2 = 0.61_wp * pt(k,j,i)
3266                         ENDIF
3267
3268                         tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
3269                                         k2 * surf_def_h(l)%qsws(m)            &
3270                                       ) * drho_air_zw(k-1)
3271                      ENDDO
3272                   ENDDO
3273!
3274!--                Treat horizontal natural surfaces
3275                   surf_s = surf_lsm_h%start_index(j,i)
3276                   surf_e = surf_lsm_h%end_index(j,i)
3277                   DO  m = surf_s, surf_e
3278                      k = surf_lsm_h%k(m)
3279
3280                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3281                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3282                         k2 = 0.61_wp * pt(k,j,i)
3283                      ELSE IF ( bulk_cloud_model )  THEN
3284                         IF ( ql(k,j,i) == 0.0_wp )  THEN
3285                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3286                            k2 = 0.61_wp * pt(k,j,i)
3287                         ELSE
3288                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3289                            temp  = theta * exner(k)
3290                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3291                                          ( q(k,j,i) - ql(k,j,i) ) *           &
3292                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3293                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3294                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3295                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3296                         ENDIF
3297                      ELSE IF ( cloud_droplets )  THEN
3298                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3299                         k2 = 0.61_wp * pt(k,j,i)
3300                      ENDIF
3301
3302                      tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
3303                                      k2 * surf_lsm_h%qsws(m)                  &
3304                                    ) * drho_air_zw(k-1)
3305                   ENDDO
3306!
3307!--                Treat horizontal urban surfaces
3308                   surf_s = surf_usm_h%start_index(j,i)
3309                   surf_e = surf_usm_h%end_index(j,i)
3310                   DO  m = surf_s, surf_e
3311                      k = surf_usm_h%k(m)
3312
3313                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3314                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3315                         k2 = 0.61_wp * pt(k,j,i)
3316                      ELSE IF ( bulk_cloud_model )  THEN
3317                         IF ( ql(k,j,i) == 0.0_wp )  THEN
3318                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3319                            k2 = 0.61_wp * pt(k,j,i)
3320                         ELSE
3321                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3322                            temp  = theta * exner(k)
3323                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3324                                          ( q(k,j,i) - ql(k,j,i) ) *           &
3325                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3326                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3327                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3328                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3329                         ENDIF
3330                      ELSE IF ( cloud_droplets )  THEN
3331                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3332                         k2 = 0.61_wp * pt(k,j,i)
3333                      ENDIF
3334
3335                      tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
3336                                      k2 * surf_usm_h%qsws(m)                  &
3337                                    ) * drho_air_zw(k-1)
3338                   ENDDO
3339
3340                ENDIF ! from IF ( use_surface_fluxes )  THEN
3341
3342                IF ( use_top_fluxes )  THEN
3343
3344                   surf_s = surf_def_h(2)%start_index(j,i)
3345                   surf_e = surf_def_h(2)%end_index(j,i)
3346                   DO  m = surf_s, surf_e
3347                      k = surf_def_h(2)%k(m)
3348
3349                      IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3350                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3351                         k2 = 0.61_wp * pt(k,j,i)
3352                      ELSE IF ( bulk_cloud_model )  THEN
3353                         IF ( ql(k,j,i) == 0.0_wp )  THEN
3354                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3355                            k2 = 0.61_wp * pt(k,j,i)
3356                         ELSE
3357                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3358                            temp  = theta * exner(k)
3359                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3360                                       ( q(k,j,i) - ql(k,j,i) ) *              &
3361                              ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3362                              ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3363                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3364                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3365                         ENDIF
3366                      ELSE IF ( cloud_droplets )  THEN
3367                         k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3368                         k2 = 0.61_wp * pt(k,j,i)
3369                      ENDIF
3370
3371                      tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
3372                                      k2 * surf_def_h(2)%qsws(m)               &
3373                                    ) * drho_air_zw(k)
3374
3375                   ENDDO
3376
3377                ENDIF ! from IF ( use_top_fluxes )  THEN
3378
3379                IF ( .NOT. diss_production )  THEN
3380
3381!--                Compute tendency for TKE-production from shear
3382                   DO  k = nzb+1, nzt
3383                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3384                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3385                                    MERGE( vpt_reference, vpt(k,j,i),          &
3386                                           use_single_reference_value ) )
3387                   ENDDO
3388
3389                ELSE
3390
3391!--                RANS mode: Compute tendency for dissipation-rate-production from shear
3392                   DO  k = nzb+1, nzt
3393                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3394                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
3395                                    MERGE( vpt_reference, vpt(k,j,i),          &
3396                                           use_single_reference_value ) ) *    &
3397                                    diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
3398                                    c_3
3399                   ENDDO
3400
3401                ENDIF
3402
3403             ENDDO
3404          ENDDO
3405
3406       ENDIF
3407
3408    ENDIF
3409
3410 END SUBROUTINE production_e
3411
3412
3413!------------------------------------------------------------------------------!
3414! Description:
3415! ------------
3416!> Production terms (shear + buoyancy) of the TKE.
3417!> Cache-optimized version
3418!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
3419!>          not considered well!
3420!------------------------------------------------------------------------------!
3421 SUBROUTINE production_e_ij( i, j, diss_production )
3422
3423    USE arrays_3d,                                                             &
3424        ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
3425
3426    USE control_parameters,                                                    &
3427        ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
3428               rho_reference, use_single_reference_value, use_surface_fluxes,  &
3429               use_top_fluxes
3430
3431    USE grid_variables,                                                        &
3432        ONLY:  ddx, dx, ddy, dy
3433
3434    USE bulk_cloud_model_mod,                                                  &
3435        ONLY:  bulk_cloud_model
3436
3437    IMPLICIT NONE
3438
3439    LOGICAL :: diss_production
3440
3441    INTEGER(iwp) ::  i       !< running index x-direction
3442    INTEGER(iwp) ::  j       !< running index y-direction
3443    INTEGER(iwp) ::  k       !< running index z-direction
3444    INTEGER(iwp) ::  l       !< running index for different surface type orientation
3445    INTEGER(iwp) ::  m       !< running index surface elements
3446    INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
3447    INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
3448    INTEGER(iwp) ::  flag_nr !< number of masking flag
3449
3450    REAL(wp)     ::  def         !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
3451    REAL(wp)     ::  flag        !< flag to mask topography
3452    REAL(wp)     ::  k1          !< temporary factor
3453    REAL(wp)     ::  k2          !< temporary factor
3454    REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
3455    REAL(wp)     ::  theta       !< virtual potential temperature
3456    REAL(wp)     ::  temp        !< theta * Exner-function
3457    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
3458    REAL(wp)     ::  usvs        !< momentum flux u"v"
3459    REAL(wp)     ::  vsus        !< momentum flux v"u"
3460    REAL(wp)     ::  wsus        !< momentum flux w"u"
3461    REAL(wp)     ::  wsvs        !< momentum flux w"v"
3462
3463    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
3464    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
3465    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
3466    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
3467    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
3468    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
3469    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
3470    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
3471    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
3472    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
3473
3474
3475
3476!
3477!-- Calculate TKE production by shear. Calculate gradients at all grid
3478!-- points first, gradients at surface-bounded grid points will be
3479!-- overwritten further below.
3480    DO  k = nzb+1, nzt
3481
3482       dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
3483       dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
3484                             u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
3485       dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
3486                             u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
3487
3488       dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
3489                             v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
3490       dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
3491       dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
3492                             v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
3493
3494       dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
3495                             w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
3496       dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
3497                             w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
3498       dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
3499
3500    ENDDO
3501
3502    flag_nr = 29
3503
3504    IF ( constant_flux_layer )  THEN
3505
3506       flag_nr = 0
3507
3508!--    Position beneath wall
3509!--    (2) - Will allways be executed.
3510!--    'bottom and wall: use u_0,v_0 and wall functions'
3511!
3512!--    Compute gradients at north- and south-facing surfaces.
3513!--    First, for default surfaces, then for urban surfaces.
3514!--    Note, so far no natural vertical surfaces implemented
3515       DO  l = 0, 1
3516          surf_s = surf_def_v(l)%start_index(j,i)
3517          surf_e = surf_def_v(l)%end_index(j,i)
3518          DO  m = surf_s, surf_e
3519             k           = surf_def_v(l)%k(m)
3520             usvs        = surf_def_v(l)%mom_flux_tke(0,m)
3521             wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
3522
3523             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
3524                             * 0.5_wp * dy
3525!
3526!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3527             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3528                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3529             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3530             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3531          ENDDO
3532!
3533!--       Natural surfaces
3534          surf_s = surf_lsm_v(l)%start_index(j,i)
3535          surf_e = surf_lsm_v(l)%end_index(j,i)
3536          DO  m = surf_s, surf_e
3537             k           = surf_lsm_v(l)%k(m)
3538             usvs        = surf_lsm_v(l)%mom_flux_tke(0,m)
3539             wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
3540
3541             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
3542                             * 0.5_wp * dy
3543!
3544!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3545             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3546                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3547             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3548             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3549          ENDDO
3550!
3551!--       Urban surfaces
3552          surf_s = surf_usm_v(l)%start_index(j,i)
3553          surf_e = surf_usm_v(l)%end_index(j,i)
3554          DO  m = surf_s, surf_e
3555             k           = surf_usm_v(l)%k(m)
3556             usvs        = surf_usm_v(l)%mom_flux_tke(0,m)
3557             wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
3558
3559             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
3560                             * 0.5_wp * dy
3561!
3562!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3563             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3564                               BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3565             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3566             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3567          ENDDO
3568       ENDDO
3569!
3570!--    Compute gradients at east- and west-facing walls
3571       DO  l = 2, 3
3572          surf_s = surf_def_v(l)%start_index(j,i)
3573          surf_e = surf_def_v(l)%end_index(j,i)
3574          DO  m = surf_s, surf_e
3575             k     = surf_def_v(l)%k(m)
3576             vsus  = surf_def_v(l)%mom_flux_tke(0,m)
3577             wsus  = surf_def_v(l)%mom_flux_tke(1,m)
3578
3579             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
3580                                * 0.5_wp * dx
3581!
3582!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3583             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3584                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3585             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3586             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3587          ENDDO
3588!
3589!--       Natural surfaces
3590          surf_s = surf_lsm_v(l)%start_index(j,i)
3591          surf_e = surf_lsm_v(l)%end_index(j,i)
3592          DO  m = surf_s, surf_e
3593             k     = surf_lsm_v(l)%k(m)
3594             vsus  = surf_lsm_v(l)%mom_flux_tke(0,m)
3595             wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
3596
3597             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
3598                                * 0.5_wp * dx
3599!
3600!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3601             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3602                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3603             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3604             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3605          ENDDO
3606!
3607!--       Urban surfaces
3608          surf_s = surf_usm_v(l)%start_index(j,i)
3609          surf_e = surf_usm_v(l)%end_index(j,i)
3610          DO  m = surf_s, surf_e
3611             k     = surf_usm_v(l)%k(m)
3612             vsus  = surf_usm_v(l)%mom_flux_tke(0,m)
3613             wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
3614
3615             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
3616                                * 0.5_wp * dx
3617!
3618!--          -1.0 for right-facing wall, 1.0 for left-facing wall
3619             sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
3620                               BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3621             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3622             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3623          ENDDO
3624       ENDDO
3625!
3626!--    Compute gradients at upward-facing surfaces
3627       surf_s = surf_def_h(0)%start_index(j,i)
3628       surf_e = surf_def_h(0)%end_index(j,i)
3629       DO  m = surf_s, surf_e
3630          k = surf_def_h(0)%k(m)
3631!
3632!--       Please note, actually, an interpolation of u_0 and v_0
3633!--       onto the grid center would be required. However, this
3634!--       would require several data transfers between 2D-grid and
3635!--       wall type. The effect of this missing interpolation is
3636!--       negligible. (See also production_e_init).
3637          dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
3638          dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
3639
3640       ENDDO
3641!
3642!--    Natural surfaces
3643       surf_s = surf_lsm_h%start_index(j,i)
3644       surf_e = surf_lsm_h%end_index(j,i)
3645       DO  m = surf_s, surf_e
3646          k = surf_lsm_h%k(m)
3647
3648          dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
3649          dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
3650
3651       ENDDO
3652!
3653!--    Urban surfaces
3654       surf_s = surf_usm_h%start_index(j,i)
3655       surf_e = surf_usm_h%end_index(j,i)
3656       DO  m = surf_s, surf_e
3657          k = surf_usm_h%k(m)
3658
3659          dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
3660          dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
3661
3662       ENDDO
3663!
3664!--    Compute gradients at downward-facing walls, only for
3665!--    non-natural default surfaces
3666       surf_s = surf_def_h(1)%start_index(j,i)
3667       surf_e = surf_def_h(1)%end_index(j,i)
3668       DO  m = surf_s, surf_e
3669          k = surf_def_h(1)%k(m)
3670
3671          dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
3672          dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
3673
3674       ENDDO
3675
3676    ENDIF
3677
3678    DO  k = nzb+1, nzt
3679
3680       def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
3681                        dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
3682                        dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
3683             2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
3684                        dwdy(k)*dvdz(k) )
3685
3686       IF ( def < 0.0_wp )  def = 0.0_wp
3687
3688       flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
3689
3690       IF ( .NOT. diss_production )  THEN
3691
3692!--       Compute tendency for TKE-production from shear
3693          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
3694
3695       ELSE
3696
3697!--       RANS mode: Compute tendency for dissipation-rate-production from shear
3698          tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
3699                        diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
3700
3701       ENDIF
3702
3703    ENDDO
3704
3705!
3706!-- If required, calculate TKE production by buoyancy
3707    IF ( .NOT. neutral )  THEN
3708
3709       IF ( .NOT. humidity )  THEN
3710
3711          IF ( ocean_mode )  THEN
3712!
3713!--          So far in the ocean no special treatment of density flux
3714!--          in the bottom and top surface layer
3715             DO  k = nzb+1, nzt
3716                tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
3717             ENDDO
3718!
3719!--          Treatment of near-surface grid points, at up- and down-
3720!--          ward facing surfaces
3721             IF ( use_surface_fluxes )  THEN
3722                DO  l = 0, 1
3723                   surf_s = surf_def_h(l)%start_index(j,i)
3724                   surf_e = surf_def_h(l)%end_index(j,i)
3725                   DO  m = surf_s, surf_e
3726                      k = surf_def_h(l)%k(m)
3727                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3728                   ENDDO
3729                ENDDO
3730             ENDIF
3731
3732             IF ( use_top_fluxes )  THEN
3733                surf_s = surf_def_h(2)%start_index(j,i)
3734                surf_e = surf_def_h(2)%end_index(j,i)
3735                DO  m = surf_s, surf_e
3736                   k = surf_def_h(2)%k(m)
3737                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3738                ENDDO
3739             ENDIF
3740
3741             IF ( .NOT. diss_production )  THEN
3742
3743!--             Compute tendency for TKE-production from shear
3744                DO  k = nzb+1, nzt
3745                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3746                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3747                                 MERGE( rho_reference, prho(k,j,i),       &
3748                                        use_single_reference_value ) )
3749                ENDDO
3750
3751             ELSE
3752
3753!--             RANS mode: Compute tendency for dissipation-rate-production from shear
3754                DO  k = nzb+1, nzt
3755                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3756                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3757                                 MERGE( rho_reference, prho(k,j,i),       &
3758                                        use_single_reference_value ) ) *  &
3759                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3760                                 c_3
3761                ENDDO
3762
3763             ENDIF
3764
3765
3766          ELSE ! or IF ( .NOT. ocean_mode )  THEN
3767
3768             DO  k = nzb+1, nzt
3769                tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
3770             ENDDO
3771
3772             IF ( use_surface_fluxes )  THEN
3773!
3774!--             Default surfaces, up- and downward-facing
3775                DO  l = 0, 1
3776                   surf_s = surf_def_h(l)%start_index(j,i)
3777                   surf_e = surf_def_h(l)%end_index(j,i)
3778                   DO  m = surf_s, surf_e
3779                      k = surf_def_h(l)%k(m)
3780                      tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3781                   ENDDO
3782                ENDDO
3783!
3784!--             Natural surfaces
3785                surf_s = surf_lsm_h%start_index(j,i)
3786                surf_e = surf_lsm_h%end_index(j,i)
3787                DO  m = surf_s, surf_e
3788                   k = surf_lsm_h%k(m)
3789                   tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
3790                ENDDO
3791!
3792!--             Urban surfaces
3793                surf_s = surf_usm_h%start_index(j,i)
3794                surf_e = surf_usm_h%end_index(j,i)
3795                DO  m = surf_s, surf_e
3796                   k = surf_usm_h%k(m)
3797                   tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
3798                ENDDO
3799             ENDIF
3800
3801             IF ( use_top_fluxes )  THEN
3802                surf_s = surf_def_h(2)%start_index(j,i)
3803                surf_e = surf_def_h(2)%end_index(j,i)
3804                DO  m = surf_s, surf_e
3805                   k = surf_def_h(2)%k(m)
3806                   tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3807                ENDDO
3808             ENDIF
3809
3810             IF ( .NOT. diss_production )  THEN
3811
3812!--             Compute tendency for TKE-production from shear
3813                DO  k = nzb+1, nzt
3814                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3815                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3816                                 MERGE( pt_reference, pt(k,j,i),          &
3817                                        use_single_reference_value ) )
3818                ENDDO
3819
3820             ELSE
3821
3822!--             RANS mode: Compute tendency for dissipation-rate-production from shear
3823                DO  k = nzb+1, nzt
3824                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3825                   tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3826                                 MERGE( pt_reference, pt(k,j,i),          &
3827                                        use_single_reference_value ) ) *  &
3828                                 diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
3829                                 c_3
3830                ENDDO
3831
3832             ENDIF
3833
3834          ENDIF ! from IF ( .NOT. ocean_mode )
3835
3836       ELSE ! or IF ( humidity )  THEN
3837
3838          DO  k = nzb+1, nzt
3839
3840             IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3841                k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3842                k2 = 0.61_wp * pt(k,j,i)
3843                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3844                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3845                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3846                                ) * dd2zu(k)
3847             ELSE IF ( bulk_cloud_model )  THEN
3848                IF ( ql(k,j,i) == 0.0_wp )  THEN
3849                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3850                   k2 = 0.61_wp * pt(k,j,i)
3851                ELSE
3852                   theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3853                   temp  = theta * exner(k)
3854                   k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
3855                                 ( q(k,j,i) - ql(k,j,i) ) *              &
3856                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
3857                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
3858                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3859                   k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3860                ENDIF
3861                tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
3862                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3863                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
3864                                ) * dd2zu(k)
3865             ELSE IF ( cloud_droplets )  THEN
3866                k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3867                k2 = 0.61_wp * pt(k,j,i)
3868                tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3869                                ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
3870                                  k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
3871                                  pt(k,j,i) * ( ql(k+1,j,i) -            &
3872                                  ql(k-1,j,i) ) ) * dd2zu(k)
3873             ENDIF
3874
3875          ENDDO
3876
3877          IF ( use_surface_fluxes )  THEN
3878!
3879!--          Treat horizontal default surfaces
3880             DO  l = 0, 1
3881                surf_s = surf_def_h(l)%start_index(j,i)
3882                surf_e = surf_def_h(l)%end_index(j,i)
3883                DO  m = surf_s, surf_e
3884                   k = surf_def_h(l)%k(m)
3885
3886                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3887                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3888                      k2 = 0.61_wp * pt(k,j,i)
3889                   ELSE IF ( bulk_cloud_model )  THEN
3890                      IF ( ql(k,j,i) == 0.0_wp )  THEN
3891                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3892                         k2 = 0.61_wp * pt(k,j,i)
3893                      ELSE
3894                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3895                         temp  = theta * exner(k)
3896                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
3897                                    ( q(k,j,i) - ql(k,j,i) ) *           &
3898                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3899                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3900                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3901                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3902                      ENDIF
3903                   ELSE IF ( cloud_droplets )  THEN
3904                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3905                      k2 = 0.61_wp * pt(k,j,i)
3906                   ENDIF
3907
3908                   tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
3909                                   k2 * surf_def_h(l)%qsws(m)            &
3910                                 ) * drho_air_zw(k-1)
3911                ENDDO
3912             ENDDO
3913!
3914!--          Treat horizontal natural surfaces
3915             surf_s = surf_lsm_h%start_index(j,i)
3916             surf_e = surf_lsm_h%end_index(j,i)
3917             DO  m = surf_s, surf_e
3918                k = surf_lsm_h%k(m)
3919
3920                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3921                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3922                   k2 = 0.61_wp * pt(k,j,i)
3923                ELSE IF ( bulk_cloud_model )  THEN
3924                   IF ( ql(k,j,i) == 0.0_wp )  THEN
3925                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3926                      k2 = 0.61_wp * pt(k,j,i)
3927                   ELSE
3928                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3929                      temp  = theta * exner(k)
3930                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3931                                    ( q(k,j,i) - ql(k,j,i) ) *           &
3932                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3933                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3934                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3935                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3936                   ENDIF
3937                ELSE IF ( cloud_droplets )  THEN
3938                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3939                   k2 = 0.61_wp * pt(k,j,i)
3940                ENDIF
3941
3942                tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
3943                                k2 * surf_lsm_h%qsws(m)                  &
3944                              ) * drho_air_zw(k-1)
3945             ENDDO
3946!
3947!--          Treat horizontal urban surfaces
3948             surf_s = surf_usm_h%start_index(j,i)
3949             surf_e = surf_usm_h%end_index(j,i)
3950             DO  m = surf_s, surf_e
3951                k = surf_usm_h%k(m)
3952
3953                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3954                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3955                   k2 = 0.61_wp * pt(k,j,i)
3956                ELSE IF ( bulk_cloud_model )  THEN
3957                   IF ( ql(k,j,i) == 0.0_wp )  THEN
3958                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3959                      k2 = 0.61_wp * pt(k,j,i)
3960                   ELSE
3961                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3962                      temp  = theta * exner(k)
3963                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
3964                                    ( q(k,j,i) - ql(k,j,i) ) *           &
3965                           ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
3966                           ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
3967                           ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3968                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3969                   ENDIF
3970                ELSE IF ( cloud_droplets )  THEN
3971                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3972                   k2 = 0.61_wp * pt(k,j,i)
3973                ENDIF
3974
3975                tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
3976                                k2 * surf_usm_h%qsws(m)                  &
3977                              ) * drho_air_zw(k-1)
3978             ENDDO
3979
3980          ENDIF ! from IF ( use_surface_fluxes )  THEN
3981
3982          IF ( use_top_fluxes )  THEN
3983
3984             surf_s = surf_def_h(2)%start_index(j,i)
3985             surf_e = surf_def_h(2)%end_index(j,i)
3986             DO  m = surf_s, surf_e
3987                k = surf_def_h(2)%k(m)
3988
3989                IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3990                   k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3991                   k2 = 0.61_wp * pt(k,j,i)
3992                ELSE IF ( bulk_cloud_model )  THEN
3993                   IF ( ql(k,j,i) == 0.0_wp )  THEN
3994                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3995                      k2 = 0.61_wp * pt(k,j,i)
3996                   ELSE
3997                      theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3998                      temp  = theta * exner(k)
3999                      k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
4000                                 ( q(k,j,i) - ql(k,j,i) ) *              &
4001                        ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
4002                        ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
4003                        ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
4004                      k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
4005                   ENDIF
4006                ELSE IF ( cloud_droplets )  THEN
4007                   k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
4008                   k2 = 0.61_wp * pt(k,j,i)
4009                ENDIF
4010
4011                tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
4012                                k2 * surf_def_h(2)%qsws(m)               &
4013                              ) * drho_air_zw(k)
4014
4015             ENDDO
4016
4017          ENDIF ! from IF ( use_top_fluxes )  THEN
4018
4019          IF ( .NOT. diss_production )  THEN
4020
4021!--          Compute tendency for TKE-production from shear
4022             DO  k = nzb+1, nzt
4023                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
4024                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
4025                              MERGE( vpt_reference, vpt(k,j,i),          &
4026                                     use_single_reference_value ) )
4027             ENDDO
4028
4029          ELSE
4030
4031!--          RANS mode: Compute tendency for dissipation-rate-production from shear
4032             DO  k = nzb+1, nzt
4033                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
4034                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
4035                              MERGE( vpt_reference, vpt(k,j,i),          &
4036                                     use_single_reference_value ) ) *    &
4037                              diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
4038                              c_3
4039             ENDDO
4040
4041          ENDIF
4042
4043       ENDIF
4044
4045    ENDIF
4046
4047 END SUBROUTINE production_e_ij
4048
4049
4050!------------------------------------------------------------------------------!
4051! Description:
4052! ------------
4053!> Diffusion and dissipation terms for the TKE.
4054!> Vector-optimized version
4055!------------------------------------------------------------------------------!
4056 SUBROUTINE diffusion_e( var, var_reference )
4057
4058    USE arrays_3d,                                                             &
4059        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
4060
4061    USE control_parameters,                                                    &
4062        ONLY:  atmos_ocean_sign, use_single_reference_value
4063
4064    USE grid_variables,                                                        &
4065        ONLY:  ddx2, ddy2
4066
4067    USE bulk_cloud_model_mod,                                                  &
4068        ONLY:  collision_turbulence
4069
4070    USE particle_attributes,                                                   &
4071        ONLY:  use_sgs_for_particles, wang_kernel
4072
4073    IMPLICIT NONE
4074
4075    INTEGER(iwp) ::  i              !< running index x direction
4076    INTEGER(iwp) ::  j              !< running index y direction
4077    INTEGER(iwp) ::  k              !< running index z direction
4078    INTEGER(iwp) ::  m              !< running index surface elements
4079
4080    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
4081    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
4082    REAL(wp)     ::  l              !< mixing length
4083    REAL(wp)     ::  var_reference  !< reference temperature
4084
4085    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
4086    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
4087
4088    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
4089
4090!
4091!-- Calculate the dissipation
4092    IF ( les_dynamic .OR. les_mw )  THEN
4093
4094       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
4095       !$ACC PRIVATE(l, l_stable, dvar_dz) &
4096       !$ACC PRESENT(diss, e, var, wall_flags_0) &
4097       !$ACC PRESENT(dd2zu, l_grid, l_wall)
4098       DO  i = nxl, nxr
4099          DO  j = nys, nyn
4100             !$ACC LOOP PRIVATE(k)
4101             DO  k = nzb+1, nzt
4102
4103                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4104                IF ( dvar_dz > 0.0_wp ) THEN
4105                   IF ( use_single_reference_value )  THEN
4106                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
4107                                 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
4108                   ELSE
4109                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
4110                                 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
4111                   ENDIF
4112                ELSE
4113                   l_stable(k) = l_grid(k)
4114                ENDIF
4115
4116             ENDDO
4117
4118             !$ACC LOOP PRIVATE(k)
4119             !DIR$ IVDEP
4120             DO  k = nzb+1, nzt
4121
4122                l  = MIN( l_wall(k,j,i), l_stable(k) )
4123
4124                diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )                &
4125                              * e(k,j,i) * SQRT( e(k,j,i) ) / l                        &
4126                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4127
4128             ENDDO
4129
4130          ENDDO
4131       ENDDO
4132
4133    ELSEIF ( rans_tke_l )  THEN
4134
4135      !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
4136      !$ACC PRIVATE(l_stable, duv2_dz2, rif, dvar_dz) &
4137      !$ACC PRESENT(diss, e, u, v, var, wall_flags_0) &
4138      !$ACC PRESENT(dd2zu, l_black, l_wall)
4139       DO  i = nxl, nxr
4140          DO  j = nys, nyn
4141!
4142!--          Calculate Richardson-flux number
4143             IF ( use_single_reference_value )  THEN
4144                !$ACC LOOP PRIVATE(k)
4145                DO  k = nzb+1, nzt
4146
4147                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4148
4149                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
4150                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
4151                              + 1E-30_wp
4152
4153                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
4154                ENDDO
4155             ELSE
4156                !$ACC LOOP PRIVATE(k)
4157                DO  k = nzb+1, nzt
4158
4159                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4160
4161                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
4162                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
4163                              + 1E-30_wp
4164
4165                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
4166                ENDDO
4167             ENDIF
4168!
4169!--          Calculate diabatic mixing length using Dyer-profile functions
4170             !$ACC LOOP PRIVATE(k)
4171             DO  k = nzb+1, nzt
4172                IF ( rif(k) >= 0.0_wp )  THEN
4173                   l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
4174                ELSE
4175                   l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
4176                ENDIF
4177             ENDDO
4178
4179             !$ACC LOOP PRIVATE(k)
4180             !DIR$ IVDEP
4181             DO  k = nzb+1, nzt
4182
4183                diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
4184                            * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4185
4186             ENDDO
4187
4188          ENDDO
4189       ENDDO
4190
4191!-- Note, in case of rans_tke_e, the dissipation is already calculated
4192!-- in prognostic_equations
4193    ENDIF
4194
4195    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
4196    !$ACC PRESENT(diss, e, km, tend, wall_flags_0) &
4197    !$ACC PRESENT(ddzu, ddzw, rho_air_zw, drho_air)
4198    DO  i = nxl, nxr
4199       DO  j = nys, nyn
4200          DO  k = nzb+1, nzt
4201
4202             tend(k,j,i) = tend(k,j,i) + (                                            &
4203                                           (                                          &
4204                       ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )            &
4205                     - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )            &
4206                                           ) * ddx2                                   &
4207                                         + (                                          &
4208                       ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )            &
4209                     - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )            &
4210                                           ) * ddy2                                   &
4211                                         + (                                          &
4212            ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)           &
4213                                                          * rho_air_zw(k)             &
4214          - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)             &
4215                                                          * rho_air_zw(k-1)           &
4216                                           ) * ddzw(k) * drho_air(k)                  &
4217                                         ) * dsig_e                                   &
4218                                           * MERGE( 1.0_wp, 0.0_wp,                   &
4219                                                    BTEST( wall_flags_0(k,j,i), 0 ) ) &
4220                          - diss(k,j,i)
4221
4222          ENDDO
4223       ENDDO
4224    ENDDO
4225
4226!
4227!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:).
4228!-- Note, bc cannot be set in tcm_boundary conditions as the dissipation
4229!-- in LES mode is only a diagnostic quantity.
4230    IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.                  &
4231         wang_kernel  .OR.  collision_turbulence  ) )  THEN
4232!
4233!--    Upward facing surfaces
4234       DO  m = 1, bc_h(0)%ns
4235          i = bc_h(0)%i(m)
4236          j = bc_h(0)%j(m)
4237          k = bc_h(0)%k(m)
4238          diss(k-1,j,i) = diss(k,j,i)
4239       ENDDO
4240!
4241!--    Downward facing surfaces
4242       DO  m = 1, bc_h(1)%ns
4243          i = bc_h(1)%i(m)
4244          j = bc_h(1)%j(m)
4245          k = bc_h(1)%k(m)
4246          diss(k+1,j,i) = diss(k,j,i)
4247       ENDDO
4248
4249    ENDIF
4250
4251 END SUBROUTINE diffusion_e
4252
4253
4254!------------------------------------------------------------------------------!
4255! Description:
4256! ------------
4257!> Diffusion and dissipation terms for the TKE.
4258!> Cache-optimized version
4259!------------------------------------------------------------------------------!
4260 SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
4261
4262    USE arrays_3d,                                                             &
4263        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
4264
4265    USE control_parameters,                                                    &
4266        ONLY:  atmos_ocean_sign, use_single_reference_value
4267
4268    USE grid_variables,                                                        &
4269        ONLY:  ddx2, ddy2
4270
4271    USE bulk_cloud_model_mod,                                                  &
4272        ONLY:  collision_turbulence
4273
4274    USE particle_attributes,                                                   &
4275        ONLY:  use_sgs_for_particles, wang_kernel
4276
4277    IMPLICIT NONE
4278
4279    INTEGER(iwp) ::  i              !< running index x direction
4280    INTEGER(iwp) ::  j              !< running index y direction
4281    INTEGER(iwp) ::  k              !< running index z direction
4282    INTEGER(iwp) ::  m              !< running index surface elements
4283    INTEGER(iwp) ::  surf_e         !< End index of surface elements at (j,i)-gridpoint
4284    INTEGER(iwp) ::  surf_s         !< Start index of surface elements at (j,i)-gridpoint
4285
4286    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
4287    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
4288    REAL(wp)     ::  l              !< mixing length
4289    REAL(wp)     ::  var_reference  !< reference temperature
4290
4291    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
4292    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
4293
4294    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
4295
4296!
4297!-- Calculate the mixing length and dissipation...
4298!-- ...in case of LES
4299    IF ( les_dynamic .OR. les_mw )  THEN
4300
4301       DO  k = nzb+1, nzt
4302          dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4303          IF ( dvar_dz > 0.0_wp ) THEN
4304             IF ( use_single_reference_value )  THEN
4305                l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
4306                            / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
4307             ELSE
4308                l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
4309                            / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
4310             ENDIF
4311          ELSE
4312             l_stable(k) = l_grid(k)
4313          ENDIF
4314       ENDDO
4315
4316       !DIR$ IVDEP
4317       DO  k = nzb+1, nzt
4318          l  = MIN( l_wall(k,j,i), l_stable(k) )
4319
4320          diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )              &
4321                        * e(k,j,i) * SQRT( e(k,j,i) ) / l                      &
4322                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4323       ENDDO
4324
4325!
4326!-- ...in case of RANS
4327    ELSEIF ( rans_tke_l )  THEN
4328
4329!
4330!--    Calculate Richardson-flux number
4331       IF ( use_single_reference_value )  THEN
4332          DO  k = nzb+1, nzt
4333             dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4334
4335             duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2           &
4336                        + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2           &
4337                        + 1E-30_wp
4338
4339             rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
4340          ENDDO
4341       ELSE
4342          DO  k = nzb+1, nzt
4343             dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4344
4345             duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2           &
4346                        + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2           &
4347                        + 1E-30_wp
4348
4349             rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp )
4350          ENDDO
4351       ENDIF
4352!
4353!--    Calculate diabatic mixing length using Dyer-profile functions
4354       DO  k = nzb+1, nzt
4355          IF ( rif(k) >= 0.0_wp )  THEN
4356             l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
4357          ELSE
4358             l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
4359          ENDIF
4360
4361       ENDDO
4362
4363       !DIR$ IVDEP
4364       DO  k = nzb+1, nzt
4365          diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
4366                      * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4367       ENDDO
4368
4369!-- Note, in case of rans_tke_e, the dissipation is already calculated
4370!-- in prognostic_equations
4371    ENDIF
4372
4373!
4374!-- Calculate the tendency term
4375    !DIR$ IVDEP
4376    DO  k = nzb+1, nzt
4377
4378       tend(k,j,i) = tend(k,j,i) + (                                           &
4379                                      (                                        &
4380                      ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )      &
4381                    - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )      &
4382                                      ) * ddx2                                 &
4383                                    + (                                        &
4384                      ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )      &
4385                    - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )      &
4386                                      ) * ddy2                                 &
4387                                    + (                                        &
4388           ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)     &
4389                                                         * rho_air_zw(k)       &
4390         - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)       &
4391                                                         * rho_air_zw(k-1)     &
4392                                      ) * ddzw(k) * drho_air(k)                &
4393                                   ) * dsig_e                                  &
4394                                     * MERGE( 1.0_wp, 0.0_wp,                  &
4395                                              BTEST( wall_flags_0(k,j,i), 0 ) )&
4396                                 - diss(k,j,i)
4397
4398    ENDDO
4399
4400!
4401!-- Set boundary conditions of dissipation if needed for calculating the sgs
4402!-- particle velocities.
4403!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
4404!-- For each surface type determine start and end index (in case of elevated
4405!-- topography several up/downward facing surfaces may exist.
4406!-- Note, bc cannot be set in tcm_boundary conditions as the dissipation
4407!-- in LES mode is only a diagnostic quantity.
4408    IF ( .NOT. rans_tke_e .AND.  ( use_sgs_for_particles  .OR.  wang_kernel    &
4409          .OR.  collision_turbulence ) )  THEN
4410       surf_s = bc_h(0)%start_index(j,i)
4411       surf_e = bc_h(0)%end_index(j,i)
4412       DO  m = surf_s, surf_e
4413          k             = bc_h(0)%k(m)
4414          diss(k-1,j,i) = diss(k,j,i)
4415       ENDDO
4416!
4417!--    Downward facing surfaces
4418       surf_s = bc_h(1)%start_index(j,i)
4419       surf_e = bc_h(1)%end_index(j,i)
4420       DO  m = surf_s, surf_e
4421          k             = bc_h(1)%k(m)
4422          diss(k+1,j,i) = diss(k,j,i)
4423       ENDDO
4424    ENDIF
4425
4426 END SUBROUTINE diffusion_e_ij
4427
4428
4429!------------------------------------------------------------------------------!
4430! Description:
4431! ------------
4432!> Diffusion term for the TKE dissipation rate
4433!> Vector-optimized version
4434!------------------------------------------------------------------------------!
4435 SUBROUTINE diffusion_diss
4436    USE arrays_3d,                                                             &
4437        ONLY:  ddzu, ddzw, drho_air, rho_air_zw
4438
4439    USE grid_variables,                                                        &
4440        ONLY:  ddx2, ddy2
4441
4442    IMPLICIT NONE
4443
4444    INTEGER(iwp) ::  i              !< running index x direction
4445    INTEGER(iwp) ::  j              !< running index y direction
4446    INTEGER(iwp) ::  k              !< running index z direction
4447
4448    REAL(wp)     ::  flag           !< flag to mask topography
4449
4450!
4451!-- Calculate the tendency terms
4452    DO  i = nxl, nxr
4453       DO  j = nys, nyn
4454          DO  k = nzb+1, nzt
4455
4456!
4457!--          Predetermine flag to mask topography
4458             flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4459
4460             tend(k,j,i) = tend(k,j,i) +                                       &
4461                         (       (                                             &
4462                 ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) )     &
4463               - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) )     &
4464                                 ) * ddx2                                      &
4465                               + (                                             &
4466                 ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) )     &
4467               - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) )     &
4468                                 ) * ddy2                                      &
4469                               + (                                             &
4470      ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1)    &
4471                                                    * rho_air_zw(k)            &
4472    - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k)      &
4473                                                    * rho_air_zw(k-1)          &
4474                                 ) * ddzw(k) * drho_air(k)                     &
4475                         ) * flag * dsig_diss                                  &
4476                         - c_2 * diss(k,j,i)**2                                &
4477                               / ( e(k,j,i) + 1.0E-20_wp ) * flag
4478
4479          ENDDO
4480       ENDDO
4481    ENDDO
4482
4483 END SUBROUTINE diffusion_diss
4484
4485
4486!------------------------------------------------------------------------------!
4487! Description:
4488! ------------
4489!> Diffusion term for the TKE dissipation rate
4490!> Cache-optimized version
4491!------------------------------------------------------------------------------!
4492 SUBROUTINE diffusion_diss_ij( i, j )
4493
4494    USE arrays_3d,                                                             &
4495        ONLY:  ddzu, ddzw, drho_air, rho_air_zw
4496
4497    USE grid_variables,                                                        &
4498        ONLY:  ddx2, ddy2
4499
4500    IMPLICIT NONE
4501
4502    INTEGER(iwp) ::  i         !< running index x direction
4503    INTEGER(iwp) ::  j         !< running index y direction
4504    INTEGER(iwp) ::  k         !< running index z direction
4505
4506    REAL(wp)     ::  flag      !< flag to mask topography
4507
4508!
4509!-- Calculate the mixing length (for dissipation)
4510    DO  k = nzb+1, nzt
4511
4512!
4513!--    Predetermine flag to mask topography
4514       flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4515
4516!
4517!--    Calculate the tendency term
4518       tend(k,j,i) =  tend(k,j,i) +                                            &
4519                   (            (                                              &
4520                ( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) )      &
4521              - ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) )      &
4522                                ) * ddx2                                       &
4523                              + (                                              &
4524                ( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) )      &
4525              - ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) )      &
4526                                ) * ddy2                                       &
4527                              + (                                              &
4528     ( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1)     &
4529                                                   * rho_air_zw(k)             &
4530   - ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k)       &
4531                                                   * rho_air_zw(k-1)           &
4532                                ) * ddzw(k) * drho_air(k)                      &
4533                   ) * flag * dsig_diss                                        &
4534                   - c_2 * diss(k,j,i)**2 / ( e(k,j,i) + 1.0E-20_wp ) * flag
4535
4536    ENDDO
4537
4538 END SUBROUTINE diffusion_diss_ij
4539
4540
4541!------------------------------------------------------------------------------!
4542! Description:
4543! ------------
4544!> Computation of the turbulent diffusion coefficients for momentum and heat.
4545!> @bug unstable stratification is not properly considered for kh in rans mode.
4546!------------------------------------------------------------------------------!
4547 SUBROUTINE tcm_diffusivities( var, var_reference )
4548
4549    USE control_parameters,                                                    &
4550        ONLY:  bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s, &
4551               e_min
4552
4553    USE surface_layer_fluxes_mod,                                              &
4554        ONLY:  phi_m
4555
4556    INTEGER(iwp) ::  i          !< loop index
4557    INTEGER(iwp) ::  j          !< loop index
4558    INTEGER(iwp) ::  k          !< loop index
4559    INTEGER(iwp) ::  m          !< loop index
4560    INTEGER(iwp) ::  n          !< loop index
4561
4562    REAL(wp) ::  var_reference  !< reference temperature
4563
4564    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
4565
4566
4567!
4568!-- Introduce an optional minimum tke
4569    IF ( e_min > 0.0_wp )  THEN
4570       !$OMP PARALLEL DO PRIVATE(i,j,k)
4571       DO  i = nxlg, nxrg
4572          DO  j = nysg, nyng
4573             DO  k = nzb+1, nzt
4574                e(k,j,i) = MAX( e(k,j,i), e_min ) *                            &
4575                        MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4576             ENDDO
4577          ENDDO
4578       ENDDO
4579    ENDIF
4580
4581!
4582!-- Call default diffusivities routine. This is always used to calculate kh.
4583    CALL tcm_diffusivities_default( var, var_reference )
4584!
4585!-- Call dynamic subgrid model to calculate km.
4586    IF ( les_dynamic )  THEN
4587       CALL tcm_diffusivities_dynamic
4588    ENDIF
4589
4590!
4591!-- In RANS mode, use MOST to calculate km and kh within the surface layer.
4592    IF ( rans_tke_e )  THEN
4593!
4594!--    Upward facing surfaces
4595!--    Default surfaces
4596       n = 0
4597       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
4598       DO  m = 1, surf_def_h(0)%ns
4599          i = surf_def_h(0)%i(m)
4600          j = surf_def_h(0)%j(m)
4601          k = surf_def_h(0)%k(m)
4602          km(k,j,i) = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) /    &
4603                      phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
4604          kh(k,j,i) = 1.35_wp * km(k,j,i)
4605       ENDDO
4606!
4607!--    Natural surfaces
4608       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
4609       DO  m = 1, surf_lsm_h%ns
4610          i = surf_lsm_h%i(m)
4611          j = surf_lsm_h%j(m)
4612          k = surf_lsm_h%k(m)
4613          km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /          &
4614                      phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
4615          kh(k,j,i) = 1.35_wp * km(k,j,i)
4616       ENDDO
4617!
4618!--    Urban surfaces
4619       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
4620       DO  m = 1, surf_usm_h%ns
4621          i = surf_usm_h%i(m)
4622          j = surf_usm_h%j(m)
4623          k = surf_usm_h%k(m)
4624          km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /          &
4625                      phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
4626          kh(k,j,i) = 1.35_wp * km(k,j,i)
4627       ENDDO
4628
4629!
4630!--    North-, south-, west and eastward facing surfaces
4631!--    Do not consider stratification at these surfaces.
4632       DO  n = 0, 3
4633!
4634!--       Default surfaces
4635          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
4636          DO  m = 1, surf_def_v(n)%ns
4637             i = surf_def_v(n)%i(m)
4638             j = surf_def_v(n)%j(m)
4639             k = surf_def_v(n)%k(m)
4640             km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m)
4641             kh(k,j,i) = 1.35_wp * km(k,j,i)
4642          ENDDO
4643!
4644!--       Natural surfaces
4645          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
4646          DO  m = 1, surf_lsm_v(n)%ns
4647             i = surf_lsm_v(n)%i(m)
4648             j = surf_lsm_v(n)%j(m)
4649             k = surf_lsm_v(n)%k(m)
4650             km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
4651             kh(k,j,i) = 1.35_wp * km(k,j,i)
4652          ENDDO
4653!
4654!--       Urban surfaces
4655          !$OMP PARALLEL DO PRIVATE(i,j,k,m)
4656          DO  m = 1, surf_usm_v(n)%ns
4657             i = surf_usm_v(n)%i(m)
4658             j = surf_usm_v(n)%j(m)
4659             k = surf_usm_v(n)%k(m)
4660             km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
4661             kh(k,j,i) = 1.35_wp * km(k,j,i)
4662          ENDDO
4663       ENDDO
4664
4665       CALL exchange_horiz( km, nbgp )
4666       CALL exchange_horiz( kh, nbgp )
4667
4668    ENDIF
4669!
4670!-- Set boundary values (Neumann conditions)
4671!-- Downward facing surfaces
4672    !$OMP PARALLEL DO PRIVATE(i,j,k)
4673    !$ACC PARALLEL LOOP PRIVATE(i,j,k) &
4674    !$ACC PRESENT(bc_h(1), kh, km)
4675    DO  m = 1, bc_h(1)%ns
4676       i = bc_h(1)%i(m)
4677       j = bc_h(1)%j(m)
4678       k = bc_h(1)%k(m)
4679       km(k+1,j,i) = km(k,j,i)
4680       kh(k+1,j,i) = kh(k,j,i)
4681    ENDDO
4682!
4683!-- Downward facing surfaces
4684    !$OMP PARALLEL DO PRIVATE(i,j,k)
4685    !$ACC PARALLEL LOOP PRIVATE(i,j,k) &
4686    !$ACC PRESENT(bc_h(0), kh, km)
4687    DO  m = 1, bc_h(0)%ns
4688       i = bc_h(0)%i(m)
4689       j = bc_h(0)%j(m)
4690       k = bc_h(0)%k(m)
4691       km(k-1,j,i) = km(k,j,i)
4692       kh(k-1,j,i) = kh(k,j,i)
4693    ENDDO
4694!
4695!-- Model top
4696    !$OMP PARALLEL DO
4697    !$ACC PARALLEL LOOP COLLAPSE(2) &
4698    !$ACC PRESENT(kh, km)
4699    DO  i = nxlg, nxrg
4700       DO  j = nysg, nyng
4701          km(nzt+1,j,i) = km(nzt,j,i)
4702          kh(nzt+1,j,i) = kh(nzt,j,i)
4703       ENDDO
4704    ENDDO
4705
4706!
4707!-- Set Neumann boundary conditions at the outflow boundaries in case of
4708!-- non-cyclic lateral boundaries
4709    IF ( bc_radiation_l )  THEN
4710       km(:,:,nxl-1) = km(:,:,nxl)
4711       kh(:,:,nxl-1) = kh(:,:,nxl)
4712    ENDIF
4713    IF ( bc_radiation_r )  THEN
4714       km(:,:,nxr+1) = km(:,:,nxr)
4715       kh(:,:,nxr+1) = kh(:,:,nxr)
4716    ENDIF
4717    IF ( bc_radiation_s )  THEN
4718       km(:,nys-1,:) = km(:,nys,:)
4719       kh(:,nys-1,:) = kh(:,nys,:)
4720    ENDIF
4721    IF ( bc_radiation_n )  THEN
4722       km(:,nyn+1,:) = km(:,nyn,:)
4723       kh(:,nyn+1,:) = kh(:,nyn,:)
4724    ENDIF
4725
4726 END SUBROUTINE tcm_diffusivities
4727
4728
4729!------------------------------------------------------------------------------!
4730! Description:
4731! ------------
4732!> Computation of the turbulent diffusion coefficients for momentum and heat
4733!> according to Prandtl-Kolmogorov.
4734!------------------------------------------------------------------------------!
4735 SUBROUTINE tcm_diffusivities_default( var, var_reference )
4736
4737    USE arrays_3d,                                                             &
4738        ONLY:  dd2zu
4739
4740    USE control_parameters,                                                    &
4741        ONLY:  atmos_ocean_sign, use_single_reference_value
4742
4743    USE statistics,                                                            &
4744        ONLY :  rmask, sums_l_l
4745
4746    IMPLICIT NONE
4747
4748    INTEGER(iwp) ::  i                   !< loop index
4749    INTEGER(iwp) ::  j                   !< loop index
4750    INTEGER(iwp) ::  k                   !< loop index
4751!$  INTEGER(iwp) ::  omp_get_thread_num  !< opemmp function to get thread number
4752    INTEGER(iwp) ::  sr                  !< statistic region
4753    INTEGER(iwp) ::  tn                  !< thread number
4754
4755    REAL(wp)     ::  duv2_dz2            !< squared vertical gradient of wind vector
4756    REAL(wp)     ::  dvar_dz             !< vertical gradient of var
4757    REAL(wp)     ::  l                   !< mixing length (single height)
4758    REAL(wp)     ::  var_reference       !< reference temperature
4759
4760    !DIR$ ATTRIBUTES ALIGN:64:: l_v, l_stable, rif
4761    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_v       !< mixing length (all heights)
4762    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
4763    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
4764
4765    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
4766
4767!
4768!-- Default thread number in case of one thread
4769    tn = 0
4770
4771!
4772!-- Initialization for calculation of the mixing length profile
4773    !$ACC KERNELS PRESENT(sums_l_l)
4774    sums_l_l = 0.0_wp
4775    !$ACC END KERNELS
4776
4777!
4778!-- Compute the turbulent diffusion coefficient for momentum
4779    !$OMP PARALLEL PRIVATE (i,j,k,l,sr,tn)
4780!$  tn = omp_get_thread_num()
4781
4782    IF ( les_dynamic .OR. les_mw )  THEN
4783       !$OMP DO
4784       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
4785       !$ACC PRIVATE(dvar_dz, l, l_stable, l_v) &
4786       !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) &
4787       !$ACC PRESENT(kh, km, sums_l_l)
4788       DO  i = nxlg, nxrg
4789          DO  j = nysg, nyng
4790             !$ACC LOOP PRIVATE(k)
4791             DO  k = nzb+1, nzt
4792!
4793!--             Determine the mixing length
4794!--             @note The following code cannot be transferred to a subroutine
4795!--             due to errors when using OpenACC directives. The execution
4796!--             crashes reliably if a subroutine is called at this point (the
4797!--             reasong for this behaviour is unknown, however).
4798                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4799                IF ( dvar_dz > 0.0_wp ) THEN
4800                   IF ( use_single_reference_value )  THEN
4801                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
4802                                  / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
4803                   ELSE
4804                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
4805                                  / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
4806                   ENDIF
4807                ELSE
4808                   l_stable(k) = l_grid(k)
4809                ENDIF
4810
4811             ENDDO
4812
4813             !$ACC LOOP PRIVATE(k)
4814             !DIR$ IVDEP
4815             DO  k = nzb+1, nzt
4816
4817                l_v(k) = MIN( l_wall(k,j,i), l_stable(k) )                      &
4818                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4819                l = l_v(k)
4820!
4821!--             Compute diffusion coefficients for momentum and heat
4822                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) )
4823                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / l_wall(k,j,i) ) * km(k,j,i)
4824
4825             ENDDO
4826!
4827!--          Summation for averaged profile (cf. flow_statistics)
4828             !$ACC LOOP PRIVATE(sr, k)
4829             DO  sr = 0, statistic_regions
4830                DO  k = nzb+1, nzt
4831                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
4832                ENDDO
4833             ENDDO
4834
4835          ENDDO
4836       ENDDO
4837
4838    ELSEIF ( rans_tke_l )  THEN
4839
4840       !$OMP DO
4841       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
4842       !$ACC PRIVATE(dvar_dz, duv2_dz2, l_stable, l_v, rif) &
4843       !$ACC PRESENT(wall_flags_0, var, dd2zu, e, u, v, l_wall, l_black, rmask) &
4844       !$ACC PRESENT(kh, km, sums_l_l)
4845       DO  i = nxlg, nxrg
4846          DO  j = nysg, nyng
4847!
4848!--          Calculate Richardson-flux number
4849             IF ( use_single_reference_value )  THEN
4850                !$ACC LOOP PRIVATE(k)
4851                !DIR$ IVDEP
4852                DO  k = nzb+1, nzt
4853                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4854
4855                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
4856                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
4857                              + 1E-30_wp
4858
4859                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
4860                ENDDO
4861             ELSE
4862                !$ACC LOOP PRIVATE(k)
4863                !DIR$ IVDEP
4864                DO  k = nzb+1, nzt
4865                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4866
4867                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
4868                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
4869                              + 1E-30_wp
4870
4871                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
4872                ENDDO
4873             ENDIF
4874!
4875!--          Calculate diabatic mixing length using Dyer-profile functions
4876!--          In case of unstable stratification, use mixing length of neutral case
4877             !$ACC LOOP PRIVATE(k)
4878             DO  k = nzb+1, nzt
4879                IF ( rif(k) >= 0.0_wp )  THEN
4880                   l_stable(k)  = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
4881                ELSE
4882                   l_stable(k)  = l_wall(k,j,i)
4883                ENDIF
4884
4885             ENDDO
4886!
4887!--          Compute diffusion coefficients for momentum and heat
4888             !$ACC LOOP PRIVATE(k)
4889             !DIR$ IVDEP
4890             DO  k = nzb+1, nzt
4891                l_v(k)    = l_stable(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4892                km(k,j,i) = c_0 * l_v(k) * SQRT( e(k,j,i) )
4893                kh(k,j,i) = km(k,j,i) / prandtl_number
4894             ENDDO
4895!
4896!--          Summation for averaged profile (cf. flow_statistics)
4897             !$ACC LOOP PRIVATE(sr, k)
4898             DO  sr = 0, statistic_regions
4899                DO  k = nzb+1, nzt
4900                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
4901                ENDDO
4902             ENDDO
4903
4904          ENDDO
4905       ENDDO
4906
4907    ELSEIF ( rans_tke_e )  THEN
4908
4909       !$OMP DO
4910       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
4911       !$ACC PRIVATE(l_v) &
4912       !$ACC PRESENT(wall_flags_0, e, diss, rmask) &
4913       !$ACC PRESENT(kh, km, sums_l_l)
4914       DO  i = nxlg, nxrg
4915          DO  j = nysg, nyng
4916!
4917!--          Compute diffusion coefficients for momentum and heat
4918             !$ACC LOOP PRIVATE(k)
4919             !DIR$ IVDEP
4920             DO  k = nzb+1, nzt
4921
4922                l_v(k) = c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / ( diss(k,j,i) + 1.0E-30_wp ) &
4923                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4924
4925                km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * l_v(k)
4926                kh(k,j,i) = km(k,j,i) / prandtl_number
4927
4928             ENDDO
4929!
4930!--          Summation for averaged profile of mixing length (cf. flow_statistics)
4931             !$ACC LOOP PRIVATE(sr, k)
4932             DO  sr = 0, statistic_regions
4933                DO  k = nzb+1, nzt
4934                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
4935                ENDDO
4936             ENDDO
4937
4938          ENDDO
4939       ENDDO
4940
4941    ENDIF
4942
4943    !$ACC KERNELS PRESENT(sums_l_l)
4944    sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn)   ! quasi boundary-condition for
4945                                                ! data output
4946    !$ACC END KERNELS
4947!$OMP END PARALLEL
4948
4949 END SUBROUTINE tcm_diffusivities_default
4950
4951
4952!------------------------------------------------------------------------------!
4953! Description:
4954! ------------
4955!> Calculates the eddy viscosity dynamically using the linear dynamic model
4956!> according to
4957!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via
4958!> stochastic analysis."
4959!> Monte Carlo Methods and Applications 14.4 (2008): 311-329.
4960!>
4961!> Furthermore dynamic bounds are used to limit the absolute value of c* as
4962!> described in
4963!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation:
4964!> Stability via realizability." Physics of Fluids 29.10 (2017): 105104.
4965!>
4966!> @author Hauke Wurps
4967!> @author Björn Maronga
4968!------------------------------------------------------------------------------!
4969 SUBROUTINE tcm_diffusivities_dynamic
4970
4971    USE arrays_3d,                                                             &
4972        ONLY:  ddzw, dzw, dd2zu, w, ug, vg
4973
4974    USE grid_variables,                                                        &
4975        ONLY : ddx, ddy, dx, dy
4976
4977    IMPLICIT NONE
4978
4979    INTEGER(iwp) ::  i           !< running index x-direction
4980    INTEGER(iwp) ::  j           !< running index y-direction
4981    INTEGER(iwp) ::  k           !< running index z-direction
4982    INTEGER(iwp) ::  l           !< running index
4983    INTEGER(iwp) ::  m           !< running index
4984
4985    REAL(wp)     ::  dudx        !< Gradient of u-component in x-direction
4986    REAL(wp)     ::  dudy        !< Gradient of u-component in y-direction
4987    REAL(wp)     ::  dudz        !< Gradient of u-component in z-direction
4988    REAL(wp)     ::  dvdx        !< Gradient of v-component in x-direction
4989    REAL(wp)     ::  dvdy        !< Gradient of v-component in y-direction
4990    REAL(wp)     ::  dvdz        !< Gradient of v-component in z-direction
4991    REAL(wp)     ::  dwdx        !< Gradient of w-component in x-direction
4992    REAL(wp)     ::  dwdy        !< Gradient of w-component in y-direction
4993    REAL(wp)     ::  dwdz        !< Gradient of w-component in z-direction
4994
4995    REAL(wp)     ::  flag        !< topography flag
4996
4997    REAL(wp)     ::  uc(-1:1,-1:1)  !< u on grid center
4998    REAL(wp)     ::  vc(-1:1,-1:1)  !< v on grid center
4999    REAL(wp)     ::  wc(-1:1,-1:1)  !< w on grid center
5000
5001    REAL(wp)     ::  ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered u
5002    REAL(wp)     ::  vt(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered v
5003    REAL(wp)     ::  wt(nzb:nzt+1,nysg:nyng,nxlg:nxrg)  !< test filtered w
5004
5005    REAL(wp)     ::  uct         !< test filtered u on grid center
5006    REAL(wp)     ::  vct         !< test filtered v on grid center
5007    REAL