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

Last change on this file since 4109 was 4109, checked in by suehring, 2 years ago

Control discretization of advection term: separate initialization of WS advection flags for momentum and scalars. In this context, resort the bits and do some minor formatting; Make initialization of scalar-advection flags more flexible, i.e. introduce an arguemnt list to indicate non-cyclic boundaries (required for decycled scalars such as chemical species or aerosols); Introduce extended 'degradation zones', where horizontal advection of passive scalars is discretized by first-order scheme at all grid points that in the vicinity of buildings (<= 3 grid points). Even though no building is within the numerical stencil, first-order scheme is used. At fourth and fifth grid point the order of the horizontal advection scheme is successively upgraded. These extended degradation zones are used to avoid stationary numerical oscillations, which are responsible for high concentration maxima that may appear under shear-free stable conditions. Therefore, an additional 3D interger array used to store control flags is introduced; Change interface for scalar advection routine; Bugfix, avoid uninitialized value sk_num in vector version of WS scalar advection; Chemistry: Decycling boundary conditions are only set at the ghost points not on the prognostic grid points; Land-surface model: Relax checks for non-consistent initialization in case static or dynamic input is provided. For example, soil_temperature or deep_soil_temperature is not mandatory any more if dynamic input is available. Also, improper settings of x_type in namelist are only checked if no static file is available.

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