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

Last change on this file since 4510 was 4510, checked in by raasch, 10 months ago

files re-formatted to follow the PALM coding standard

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