Home

Context Navigation

turbulence_closure_mod.f90 @ 3121

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