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

Last change on this file since 4486 was 4486, checked in by maronga, 18 months ago

1.5-order-dai SGS schema modified for use with topography

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