Home

Context Navigation

turbulence_closure_mod.f90 @ 3083

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