Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |