Home

Context Navigation

turbulence_closure_mod.f90 @ 3430

Last change on this file since 3430 was 3430, checked in by maronga, 6 years ago
added support for buildings with dynamic sgs model
Property svn:keywords set to `Id`
File size: 197.0 KB

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