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

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

last changes documented

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