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

Last change on this file since 4180 was 4180, checked in by scharf, 22 months ago

removed comments in 'Former revisions' section that are older than 01.01.2019

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