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

Last change on this file since 4370 was 4370, checked in by raasch, 4 years ago

bugfixes for previous commit: unused variables removed, Temperton-fft usage on GPU, openacc porting of vector version of Obukhov length calculation, collective read switched off on NEC to avoid hanging; some vector directives added in prognostic equations to force vectorization on Intel19 compiler, configuration files for NEC Aurora added

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