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

Last change on this file since 4168 was 4168, checked in by suehring, 23 months ago

Replace get_topography_top_index functions by pre-calculated arrays in order to save computational resources

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