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

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