Home

Context Navigation

turbulence_closure_mod.f90 @ 3120

Last change on this file since 3120 was 3120, checked in by gronemeier, 5 years ago
implementation of dynamic sgs model
Property svn:keywords set to `Id`
File size: 219.2 KB

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