Home

Context Navigation

turbulence_closure_mod.f90 @ 3129

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