Home

Context Navigation

turbulence_closure_mod.f90 @ 4495

Last change on this file since 4495 was 4495, checked in by raasch, 5 years ago
restart data handling with MPI-IO added, first part
Property svn:keywords set to `Id`
File size: 215.8 KB

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