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

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

Bugfix concerning ACC directive for non-allocated array in turbulence_closure_mod; test case results updated

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