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

Last change on this file since 4586 was 4586, checked in by gronemeier, 3 months ago

renamed Richardson flux number into gradient Richardson number (model_1d_mod.f90, turbulence_closure_mod.f90, header.f90, surface_mod.f90) and zeta (header.f90);
do not add whitespaces at current-revision section (document_changes)

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