Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |