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

Last change on this file since 4495 was 4495, checked in by raasch, 18 months ago

restart data handling with MPI-IO added, first part

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