Home

Context Navigation

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

Last change on this file since 4771 was 4768, checked in by suehring, 5 years ago
Enable 3D data output also with 64-bit precision
Property svn:keywords set to `Id`
File size: 234.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |