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

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

add comment

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