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

Last change on this file since 4481 was 4481, checked in by maronga, 18 months ago

Bugfix for copyright updates in document_changes; copyright update applied to all files

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