Home

Context Navigation

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

Last change on this file since 4102 was 4102, checked in by suehring, 6 years ago
Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.
Property svn:keywords set to `Id`
File size: 204.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |