source: palm/trunk/SOURCE/turbulence_closure_mod.f90

Last change on this file was 4835, checked in by raasch, 6 months ago

openmp bugfix (some private statements were missing); output format for cpu measures slightly changed

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