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

Last change on this file since 4346 was 4346, checked in by motisi, 23 months ago

Introduction of wall_flags_total_0, which currently sets bits based on static topography information used in wall_flags_static_0

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