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

Last change on this file since 4102 was 4102, checked in by suehring, 2 years ago

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

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