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

Last change on this file since 4473 was 4473, checked in by gronemeier, 3 years ago

revised naming of mixing length and some further cleaning of mixing-length calculation

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