Home

Context Navigation

source: palm/trunk/SOURCE/turbulence_closure_mod.f90 @ 3398

Last change on this file since 3398 was 3398, checked in by knoop, 6 years ago
Refactored production_e and production_e_ij (removed code redundancy)
Property svn:keywords set to `Id`
File size: 196.7 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 3398 2018-10-22 19:30:24Z knoop $
27	! Refactored production_e and production_e_ij (removed code redundancy)
28	!
29	! 3386 2018-10-19 16:28:22Z gronemeier
30	! Renamed tcm_prognostic to tcm_prognostic_equations
31	!
32	! 3385 2018-10-19 14:52:29Z knoop
33	! Restructured loops and ifs in production_e to ease vectorization on GPUs
34	!
35	! 3300 2018-10-02 14:16:54Z gronemeier
36	! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall
37	! calculation
38	! - removed maxloc call as this produced different results for different
39	! compiler options
40	!
41	! 3294 2018-10-01 02:37:10Z raasch
42	! changes concerning modularization of ocean option
43	!
44	! 3274 2018-09-24 15:42:55Z knoop
45	! Modularization of all bulk cloud physics code components
46	!
47	! 3245 2018-09-13 14:08:16Z knoop
48	! unused variables removed, shortest_distance has wp now
49	!
50	! 3183 2018-07-27 14:25:55Z suehring
51	! Rename variables and remove unused variable from USE statement
52	!
53	! 3182 2018-07-27 13:36:03Z suehring
54	! Use MOST for km only in RANS mode
55	!
56	! 3130 2018-07-16 11:08:55Z gronemeier
57	! - move boundary condition of km and kh to tcm_diffusivities
58	! - calculate km at boundaries according to MOST
59	! - move phi_m to surface_layer_fluxes_mod
60	!
61	! 3129 2018-07-16 07:45:13Z gronemeier
62	! - move limitation of diss to boundary_conds
63	! - move boundary conditions for e and diss to boundary_conds
64	! - consider non-default surfaces in tcm_diffusivities
65	! - use z_mo within surface layer instead of calculating it
66	! - resort output after case select -> reduced code duplication
67	! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d
68	!
69	! 3121 2018-07-11 18:46:49Z gronemeier
70	! - created the function phi_m
71	! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
72	! surfaces
73	!
74	! 3120 2018-07-11 18:30:57Z gronemeier
75	! - changed tcm_diffusivities to tcm_diffusivities_default
76	! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
77	! and tcm_diffusivities_dynamic
78	!
79	! 3086 2018-06-25 09:08:04Z gronemeier
80	! bugfix: set rans_const_sigma(1) = 1.3
81	!
82	! 3083 2018-06-19 14:03:12Z gronemeier
83	! - set limits of diss at the end of prognostic equations
84	! - call production_e to calculate production term of diss
85	! - limit change of diss to -90% to +100%
86	! - remove factor 0.5 from diffusion_diss_ij
87	! - rename c_m into c_0, and c_h into c_4
88	! - add rans_const_c and rans_const_sigma as namelist parameters
89	! - add calculation of mixing length for profile output in case of rans_tke_e
90	! - changed format of annotations to comply with doxygen standards
91	! - calculate and save dissipation rate during rans_tke_l mode
92	! - set bc at vertical walls for e, diss, km, kh
93	! - bugfix: set l_wall = 0.0 within buildings
94	! - set l_wall at bottom and top boundary (rans-mode)
95	! - bugfix in production term for dissipation rate
96	! - bugfix in diffusion of dissipation rate
97	! - disable check for 1D model if rans_tke_e is used
98	! - bugfixes for initialization (rans-mode):
99	! - correction of dissipation-rate formula
100	! - calculate km based on l_wall
101	! - initialize diss if 1D model is not used
102	!
103	! 3045 2018-05-28 07:55:41Z Giersch
104	! Error message revised
105	!
106	! 3014 2018-05-09 08:42:38Z maronga
107	! Bugfix: nzb_do and nzt_do were not used for 3d data output
108	!
109	! 3004 2018-04-27 12:33:25Z Giersch
110	! Further allocation checks implemented
111	!
112	! 2938 2018-03-27 15:52:42Z suehring
113	! Further todo's
114	!
115	! 2936 2018-03-27 14:49:27Z gronemeier
116	! - defined l_grid only within this module
117	! - Moved l_wall definition from modules.f90
118	! - Get level of highest topography, used to limit upward distance calculation
119	! - Consider cyclic boundary conditions for mixing length calculation
120	! - Moved copy of wall_flags into subarray to subroutine
121	! - Implemented l_wall calculation in case of RANS simulation
122	! - Moved init of l_black to tcm_init_mixing_length
123	! - Moved init_mixing_length from init_grid.f90 and
124	! renamed it to tcm_init_mixing_length
125	!
126	! 2764 2018-01-22 09:25:36Z gronemeier
127	! Bugfix: remove duplicate SAVE statements
128	!
129	! 2746 2018-01-15 12:06:04Z suehring
130	! Move flag plant canopy to modules
131	!
132	! 2718 2018-01-02 08:49:38Z maronga
133	! Corrected "Former revisions" section
134	!
135	! 2701 2017-12-15 15:40:50Z suehring
136	! Changes from last commit documented
137	!
138	! 2698 2017-12-14 18:46:24Z suehring
139	! Bugfix in get_topography_top_index
140	!
141	! 2696 2017-12-14 17:12:51Z kanani
142	! Initial revision
143	!
144	!
145	!
146	!
147	! Authors:
148	! --------
149	! @author Tobias Gronemeier
150	! @author Hauke Wurps
151	!
152	! Description:
153	! ------------
154	!> This module contains the available turbulence closures for PALM.
155	!>
156	!>
157	!> @todo test initialization for all possibilities
158	!> add OpenMP directives whereever possible
159	!> remove debug output variables (dummy1, dummy2, dummy3)
160	!> @todo Check for random disturbances
161	!> @note <Enter notes on the module>
162	!------------------------------------------------------------------------------!
163	MODULE turbulence_closure_mod
164
165
166	#if defined( __nopointer )
167	USE arrays_3d, &
168	ONLY: diss, diss_p, dzu, e, e_p, kh, km, &
169	mean_inflow_profiles, prho, pt, tdiss_m, te_m, tend, u, v, vpt, w
170	#else
171	USE arrays_3d, &
172	ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, &
173	e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, &
174	te_m, tend, u, v, vpt, w
175	#endif
176
177	USE basic_constants_and_equations_mod, &
178	ONLY: g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
179
180	USE control_parameters, &
181	ONLY: constant_diffusion, dt_3d, e_init, humidity, &
182	initializing_actions, intermediate_timestep_count, &
183	intermediate_timestep_count_max, km_constant, &
184	les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
185	pt_reference, rans_mode, rans_tke_e, rans_tke_l, &
186	simulated_time,timestep_scheme, turbulence_closure, &
187	turbulent_inflow, use_upstream_for_tke, vpt_reference, &
188	ws_scheme_sca
189
190	USE advec_ws, &
191	ONLY: advec_s_ws
192
193	USE advec_s_bc_mod, &
194	ONLY: advec_s_bc
195
196	USE advec_s_pw_mod, &
197	ONLY: advec_s_pw
198
199	USE advec_s_up_mod, &
200	ONLY: advec_s_up
201
202	USE cpulog, &
203	ONLY: cpu_log, log_point, log_point_s
204
205	USE indices, &
206	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
207	wall_flags_0
208
209	USE kinds
210
211	USE ocean_mod, &
212	ONLY: prho_reference
213
214	USE pegrid
215
216	USE plant_canopy_model_mod, &
217	ONLY: pcm_tendency
218
219	USE statistics, &
220	ONLY: hom, hom_sum, statistic_regions
221
222	USE user_actions_mod, &
223	ONLY: user_actions
224
225
226	IMPLICIT NONE
227
228
229	REAL(wp) :: c_0 !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
230	REAL(wp) :: c_1 !< model constant for RANS mode
231	REAL(wp) :: c_2 !< model constant for RANS mode
232	REAL(wp) :: c_3 !< model constant for RANS mode
233	REAL(wp) :: c_4 !< model constant for RANS mode
234	REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length
235	REAL(wp) :: dsig_e = 1.0_wp !< factor to calculate Ke from Km (1/sigma_e)
236	REAL(wp) :: dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
237
238	REAL(wp), DIMENSION(0:4) :: rans_const_c = & !< model constants for RANS mode (namelist param)
239	(/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
240
241	REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
242	(/ 1.0_wp, 1.30_wp /)
243
244	REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length according to Blackadar
245
246	REAL(wp), DIMENSION(:), ALLOCATABLE :: l_grid !< geometric mean of grid sizes dx, dy, dz
247
248	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< near-wall mixing length
249
250	!> @todo remove debug variables
251	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
252	diss_prod1, diss_adve1, diss_diff1, &
253	diss_prod2, diss_adve2, diss_diff2, &
254	diss_prod3, diss_adve3, diss_diff3, &
255	dummy1, dummy2, dummy3
256
257
258	PUBLIC c_0, rans_const_c, rans_const_sigma
259
260	!
261	!-- PALM interfaces:
262	!-- Input parameter checks to be done in check_parameters
263	INTERFACE tcm_check_parameters
264	MODULE PROCEDURE tcm_check_parameters
265	END INTERFACE tcm_check_parameters
266
267	!
268	!-- Data output checks for 2D/3D data to be done in check_parameters
269	INTERFACE tcm_check_data_output
270	MODULE PROCEDURE tcm_check_data_output
271	END INTERFACE tcm_check_data_output
272
273	!
274	!-- Definition of data output quantities
275	INTERFACE tcm_define_netcdf_grid
276	MODULE PROCEDURE tcm_define_netcdf_grid
277	END INTERFACE tcm_define_netcdf_grid
278
279	!
280	!-- Averaging of 3D data for output
281	INTERFACE tcm_3d_data_averaging
282	MODULE PROCEDURE tcm_3d_data_averaging
283	END INTERFACE tcm_3d_data_averaging
284
285	!
286	!-- Data output of 2D quantities
287	INTERFACE tcm_data_output_2d
288	MODULE PROCEDURE tcm_data_output_2d
289	END INTERFACE tcm_data_output_2d
290
291	!
292	!-- Data output of 3D data
293	INTERFACE tcm_data_output_3d
294	MODULE PROCEDURE tcm_data_output_3d
295	END INTERFACE tcm_data_output_3d
296
297	!
298	!-- Initialization actions
299	INTERFACE tcm_init
300	MODULE PROCEDURE tcm_init
301	END INTERFACE tcm_init
302
303	!
304	!-- Initialization of arrays
305	INTERFACE tcm_init_arrays
306	MODULE PROCEDURE tcm_init_arrays
307	END INTERFACE tcm_init_arrays
308
309	!
310	!-- Initialization of TKE production term
311	INTERFACE production_e_init
312	MODULE PROCEDURE production_e_init
313	END INTERFACE production_e_init
314
315	!
316	!-- Prognostic equations for TKE and TKE dissipation rate
317	INTERFACE tcm_prognostic_equations
318	MODULE PROCEDURE tcm_prognostic_equations
319	MODULE PROCEDURE tcm_prognostic_equations_ij
320	END INTERFACE tcm_prognostic_equations
321
322	!
323	!-- Production term for TKE
324	INTERFACE production_e
325	MODULE PROCEDURE production_e
326	MODULE PROCEDURE production_e_ij
327	END INTERFACE production_e
328
329	!
330	!-- Diffusion term for TKE
331	INTERFACE diffusion_e
332	MODULE PROCEDURE diffusion_e
333	MODULE PROCEDURE diffusion_e_ij
334	END INTERFACE diffusion_e
335
336	!
337	!-- Diffusion term for TKE dissipation rate
338	INTERFACE diffusion_diss
339	MODULE PROCEDURE diffusion_diss
340	MODULE PROCEDURE diffusion_diss_ij
341	END INTERFACE diffusion_diss
342
343	!
344	!-- Mixing length for LES case
345	INTERFACE mixing_length_les
346	MODULE PROCEDURE mixing_length_les
347	END INTERFACE mixing_length_les
348
349	!
350	!-- Mixing length for RANS case
351	INTERFACE mixing_length_rans
352	MODULE PROCEDURE mixing_length_rans
353	END INTERFACE mixing_length_rans
354
355	!
356	!-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic
357	INTERFACE tcm_diffusivities
358	MODULE PROCEDURE tcm_diffusivities
359	END INTERFACE tcm_diffusivities
360
361	!
362	!-- Calculate diffusivities
363	INTERFACE tcm_diffusivities_default
364	MODULE PROCEDURE tcm_diffusivities_default
365	END INTERFACE tcm_diffusivities_default
366
367	!
368	!-- Calculate diffusivities according to dynamic sgs model
369	INTERFACE tcm_diffusivities_dynamic
370	MODULE PROCEDURE tcm_diffusivities_dynamic
371	END INTERFACE tcm_diffusivities_dynamic
372
373	!
374	!-- Box-filter method for dynamic sgs model
375	INTERFACE tcm_box_filter_2d
376	MODULE PROCEDURE tcm_box_filter_2d_single
377	MODULE PROCEDURE tcm_box_filter_2d_array
378	END INTERFACE tcm_box_filter_2d
379
380	!
381	!-- Swapping of time levels (required for prognostic variables)
382	INTERFACE tcm_swap_timelevel
383	MODULE PROCEDURE tcm_swap_timelevel
384	END INTERFACE tcm_swap_timelevel
385
386	SAVE
387
388	PRIVATE
389	!
390	!-- Add INTERFACES that must be available to other modules (alphabetical order)
391	PUBLIC production_e_init, tcm_3d_data_averaging, tcm_check_data_output, &
392	tcm_check_parameters, tcm_data_output_2d, tcm_data_output_3d, &
393	tcm_define_netcdf_grid, tcm_diffusivities, tcm_init, &
394	tcm_init_arrays, tcm_prognostic_equations, tcm_swap_timelevel
395
396
397	CONTAINS
398
399	!------------------------------------------------------------------------------!
400	! Description:
401	! ------------
402	!> Check parameters routine for turbulence closure module.
403	!> @todo remove rans_mode from initialization namelist and rework checks
404	!> The way it is implemented at the moment, the user has to set two variables
405	!> so that the RANS mode is working. It would be better if only one parameter
406	!> has to be set.
407	!> 2018-06-18, gronemeier
408	!------------------------------------------------------------------------------!
409	SUBROUTINE tcm_check_parameters
410
411	USE control_parameters, &
412	ONLY: message_string, turbulent_inflow, turbulent_outflow
413
414	IMPLICIT NONE
415
416	!
417	!-- Define which turbulence closure is going to be used
418	IF ( rans_mode ) THEN
419
420	!
421	!-- Assign values to constants for RANS mode
422	dsig_e = 1.0_wp / rans_const_sigma(1)
423	dsig_diss = 1.0_wp / rans_const_sigma(2)
424
425	c_0 = rans_const_c(0)
426	c_1 = rans_const_c(1)
427	c_2 = rans_const_c(2)
428	c_3 = rans_const_c(3) !> @todo clarify how to switch between different models
429	c_4 = rans_const_c(4)
430
431	SELECT CASE ( TRIM( turbulence_closure ) )
432
433	CASE ( 'TKE-l' )
434	rans_tke_l = .TRUE.
435
436	CASE ( 'TKE-e' )
437	rans_tke_e = .TRUE.
438
439	CASE DEFAULT
440	message_string = 'Unknown turbulence closure: ' // &
441	TRIM( turbulence_closure )
442	CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
443
444	END SELECT
445
446	IF ( turbulent_inflow .OR. turbulent_outflow ) THEN
447	message_string = 'turbulent inflow/outflow is not yet '// &
448	'implemented for RANS mode'
449	CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
450	ENDIF
451
452	message_string = 'RANS mode is still in development! ' // &
453	'&Not all features of PALM are yet compatible '// &
454	'with RANS mode. &Use at own risk!'
455	CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 )
456
457	ELSE
458
459	c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980)
460
461	dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
462	!of K_e which is used in RANS mode
463
464	SELECT CASE ( TRIM( turbulence_closure ) )
465
466	CASE ( 'Moeng_Wyngaard' )
467	les_mw = .TRUE.
468
469	CASE ( 'dynamic' )
470	les_dynamic = .TRUE.
471
472	CASE DEFAULT
473	!> @todo rework this part so that only one call of this error exists
474	message_string = 'Unknown turbulence closure: ' // &
475	TRIM( turbulence_closure )
476	CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
477
478	END SELECT
479
480	ENDIF
481
482	END SUBROUTINE tcm_check_parameters
483
484	!------------------------------------------------------------------------------!
485	! Description:
486	! ------------
487	!> Check data output.
488	!------------------------------------------------------------------------------!
489	SUBROUTINE tcm_check_data_output( var, unit )
490
491	IMPLICIT NONE
492
493	CHARACTER (LEN=*) :: unit !< unit of output variable
494	CHARACTER (LEN=*) :: var !< name of output variable
495
496
497	SELECT CASE ( TRIM( var ) )
498
499	CASE ( 'diss' )
500	unit = 'm2/s3'
501
502	CASE ( 'diss1', 'diss2', & !> @todo remove later
503	'diss_prod1', 'diss_adve1', 'diss_diff1', &
504	'diss_prod2', 'diss_adve2', 'diss_diff2', &
505	'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3' )
506	unit = 'debug output'
507
508	CASE ( 'kh', 'km' )
509	unit = 'm2/s'
510
511	CASE DEFAULT
512	unit = 'illegal'
513
514	END SELECT
515
516	END SUBROUTINE tcm_check_data_output
517
518
519	!------------------------------------------------------------------------------!
520	! Description:
521	! ------------
522	!> Define appropriate grid for netcdf variables.
523	!> It is called out from subroutine netcdf.
524	!------------------------------------------------------------------------------!
525	SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
526
527	IMPLICIT NONE
528
529	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable
530	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable
531	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable
532	CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable
533
534	LOGICAL, INTENT(OUT) :: found !< flag if output variable is found
535
536	found = .TRUE.
537
538	!
539	!-- Check for the grid
540	SELECT CASE ( TRIM( var ) )
541
542	CASE ( 'diss', 'diss_xy', 'diss_xz', 'diss_yz' )
543	grid_x = 'x'
544	grid_y = 'y'
545	grid_z = 'zu'
546
547	CASE ( 'diss1', 'diss2', & !> @todo remove later
548	'diss_prod1', 'diss_adve1', 'diss_diff1', &
549	'diss_prod2', 'diss_adve2', 'diss_diff2', &
550	'diss_prod3', 'diss_adve3', 'diss_diff3', 'dummy3' )
551	grid_x = 'x'
552	grid_y = 'y'
553	grid_z = 'zu'
554
555	CASE ( 'kh', 'kh_xy', 'kh_xz', 'kh_yz' )
556	grid_x = 'x'
557	grid_y = 'y'
558	grid_z = 'zu'
559
560	CASE ( 'km', 'km_xy', 'km_xz', 'km_yz' )
561	grid_x = 'x'
562	grid_y = 'y'
563	grid_z = 'zu'
564
565	CASE DEFAULT
566	found = .FALSE.
567	grid_x = 'none'
568	grid_y = 'none'
569	grid_z = 'none'
570
571	END SELECT
572
573	END SUBROUTINE tcm_define_netcdf_grid
574
575
576	!------------------------------------------------------------------------------!
577	! Description:
578	! ------------
579	!> Average 3D data.
580	!------------------------------------------------------------------------------!
581	SUBROUTINE tcm_3d_data_averaging( mode, variable )
582
583
584	USE averaging, &
585	ONLY: diss_av, kh_av, km_av
586
587	USE control_parameters, &
588	ONLY: average_count_3d
589
590	IMPLICIT NONE
591
592	CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average'
593	CHARACTER (LEN=*) :: variable !< name of variable
594
595	INTEGER(iwp) :: i !< loop index
596	INTEGER(iwp) :: j !< loop index
597	INTEGER(iwp) :: k !< loop index
598
599	IF ( mode == 'allocate' ) THEN
600
601	SELECT CASE ( TRIM( variable ) )
602
603	CASE ( 'diss' )
604	IF ( .NOT. ALLOCATED( diss_av ) ) THEN
605	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
606	ENDIF
607	diss_av = 0.0_wp
608
609	CASE ( 'kh' )
610	IF ( .NOT. ALLOCATED( kh_av ) ) THEN
611	ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
612	ENDIF
613	kh_av = 0.0_wp
614
615	CASE ( 'km' )
616	IF ( .NOT. ALLOCATED( km_av ) ) THEN
617	ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
618	ENDIF
619	km_av = 0.0_wp
620
621	CASE DEFAULT
622	CONTINUE
623
624	END SELECT
625
626	ELSEIF ( mode == 'sum' ) THEN
627
628	SELECT CASE ( TRIM( variable ) )
629
630	CASE ( 'diss' )
631	IF ( ALLOCATED( diss_av ) ) THEN
632	DO i = nxlg, nxrg
633	DO j = nysg, nyng
634	DO k = nzb, nzt+1
635	diss_av(k,j,i) = diss_av(k,j,i) + diss(k,j,i)
636	ENDDO
637	ENDDO
638	ENDDO
639	ENDIF
640
641	CASE ( 'kh' )
642	IF ( ALLOCATED( kh_av ) ) THEN
643	DO i = nxlg, nxrg
644	DO j = nysg, nyng
645	DO k = nzb, nzt+1
646	kh_av(k,j,i) = kh_av(k,j,i) + kh(k,j,i)
647	ENDDO
648	ENDDO
649	ENDDO
650	ENDIF
651
652	CASE ( 'km' )
653	IF ( ALLOCATED( km_av ) ) THEN
654	DO i = nxlg, nxrg
655	DO j = nysg, nyng
656	DO k = nzb, nzt+1
657	km_av(k,j,i) = km_av(k,j,i) + km(k,j,i)
658	ENDDO
659	ENDDO
660	ENDDO
661	ENDIF
662
663	CASE DEFAULT
664	CONTINUE
665
666	END SELECT
667
668	ELSEIF ( mode == 'average' ) THEN
669
670	SELECT CASE ( TRIM( variable ) )
671
672	CASE ( 'diss' )
673	IF ( ALLOCATED( diss_av ) ) THEN
674	DO i = nxlg, nxrg
675	DO j = nysg, nyng
676	DO k = nzb, nzt+1
677	diss_av(k,j,i) = diss_av(k,j,i) &
678	/ REAL( average_count_3d, KIND=wp )
679	ENDDO
680	ENDDO
681	ENDDO
682	ENDIF
683
684	CASE ( 'kh' )
685	IF ( ALLOCATED( kh_av ) ) THEN
686	DO i = nxlg, nxrg
687	DO j = nysg, nyng
688	DO k = nzb, nzt+1
689	kh_av(k,j,i) = kh_av(k,j,i) &
690	/ REAL( average_count_3d, KIND=wp )
691	ENDDO
692	ENDDO
693	ENDDO
694	ENDIF
695
696	CASE ( 'km' )
697	IF ( ALLOCATED( km_av ) ) THEN
698	DO i = nxlg, nxrg
699	DO j = nysg, nyng
700	DO k = nzb, nzt+1
701	km_av(k,j,i) = km_av(k,j,i) &
702	/ REAL( average_count_3d, KIND=wp )
703	ENDDO
704	ENDDO
705	ENDDO
706	ENDIF
707
708	END SELECT
709
710	ENDIF
711
712	END SUBROUTINE tcm_3d_data_averaging
713
714
715	!------------------------------------------------------------------------------!
716	! Description:
717	! ------------
718	!> Define 2D output variables.
719	!------------------------------------------------------------------------------!
720	SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, &
721	nzb_do, nzt_do )
722
723	USE averaging, &
724	ONLY: diss_av, kh_av, km_av
725
726	IMPLICIT NONE
727
728	CHARACTER (LEN=*) :: grid !< name of vertical grid
729	CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz'
730	CHARACTER (LEN=*) :: variable !< name of variable
731
732	INTEGER(iwp) :: av !< flag for (non-)average output
733	INTEGER(iwp) :: flag_nr !< number of masking flag
734	INTEGER(iwp) :: i !< loop index
735	INTEGER(iwp) :: j !< loop index
736	INTEGER(iwp) :: k !< loop index
737	INTEGER(iwp) :: nzb_do !< vertical output index (bottom)
738	INTEGER(iwp) :: nzt_do !< vertical output index (top)
739
740	LOGICAL :: found !< flag if output variable is found
741	LOGICAL :: resorted !< flag if output is already resorted
742
743	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
744
745	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
746	!< array to which output data is resorted to
747
748	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
749
750	found = .TRUE.
751	resorted = .FALSE.
752	!
753	!-- Set masking flag for topography for not resorted arrays
754	flag_nr = 0
755
756	SELECT CASE ( TRIM( variable ) )
757
758	CASE ( 'diss_xy', 'diss_xz', 'diss_yz' )
759	IF ( av == 0 ) THEN
760	to_be_resorted => diss
761	ELSE
762	IF ( .NOT. ALLOCATED( diss_av ) ) THEN
763	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
764	diss_av = REAL( fill_value, KIND = wp )
765	ENDIF
766	to_be_resorted => diss_av
767	ENDIF
768	IF ( mode == 'xy' ) grid = 'zu'
769
770	CASE ( 'kh_xy', 'kh_xz', 'kh_yz' )
771	IF ( av == 0 ) THEN
772	to_be_resorted => kh
773	ELSE
774	IF ( .NOT. ALLOCATED( kh_av ) ) THEN
775	ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
776	kh_av = REAL( fill_value, KIND = wp )
777	ENDIF
778	to_be_resorted => kh_av
779	ENDIF
780	IF ( mode == 'xy' ) grid = 'zu'
781
782	CASE ( 'km_xy', 'km_xz', 'km_yz' )
783	IF ( av == 0 ) THEN
784	to_be_resorted => km
785	ELSE
786	IF ( .NOT. ALLOCATED( km_av ) ) THEN
787	ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
788	km_av = REAL( fill_value, KIND = wp )
789	ENDIF
790	to_be_resorted => km_av
791	ENDIF
792	IF ( mode == 'xy' ) grid = 'zu'
793
794	CASE DEFAULT
795	found = .FALSE.
796	grid = 'none'
797
798	END SELECT
799
800	IF ( found .AND. .NOT. resorted ) THEN
801	DO i = nxl, nxr
802	DO j = nys, nyn
803	DO k = nzb_do, nzt_do
804	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
805	REAL( fill_value, KIND = wp ), &
806	BTEST( wall_flags_0(k,j,i), flag_nr ) )
807	ENDDO
808	ENDDO
809	ENDDO
810	ENDIF
811
812	END SUBROUTINE tcm_data_output_2d
813
814
815	!------------------------------------------------------------------------------!
816	! Description:
817	! ------------
818	!> Define 3D output variables.
819	!------------------------------------------------------------------------------!
820	SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
821
822
823	USE averaging, &
824	ONLY: diss_av, kh_av, km_av
825
826	IMPLICIT NONE
827
828	CHARACTER (LEN=*) :: variable !< name of variable
829
830	INTEGER(iwp) :: av !< flag for (non-)average output
831	INTEGER(iwp) :: flag_nr !< number of masking flag
832	INTEGER(iwp) :: i !< loop index
833	INTEGER(iwp) :: j !< loop index
834	INTEGER(iwp) :: k !< loop index
835	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
836	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
837
838	LOGICAL :: found !< flag if output variable is found
839	LOGICAL :: resorted !< flag if output is already resorted
840
841	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
842
843	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
844	!< array to which output data is resorted to
845
846	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
847
848	found = .TRUE.
849	resorted = .FALSE.
850	!
851	!-- Set masking flag for topography for not resorted arrays
852	flag_nr = 0
853
854	SELECT CASE ( TRIM( variable ) )
855
856	CASE ( 'diss' )
857	IF ( av == 0 ) THEN
858	to_be_resorted => diss
859	ELSE
860	IF ( .NOT. ALLOCATED( diss_av ) ) THEN
861	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
862	diss_av = REAL( fill_value, KIND = wp )
863	ENDIF
864	to_be_resorted => diss_av
865	ENDIF
866
867	CASE ( 'kh' )
868	IF ( av == 0 ) THEN
869	to_be_resorted => kh
870	ELSE
871	IF ( .NOT. ALLOCATED( kh_av ) ) THEN
872	ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
873	kh_av = REAL( fill_value, KIND = wp )
874	ENDIF
875	to_be_resorted => kh_av
876	ENDIF
877
878	CASE ( 'km' )
879	IF ( av == 0 ) THEN
880	to_be_resorted => km
881	ELSE
882	IF ( .NOT. ALLOCATED( km_av ) ) THEN
883	ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
884	km_av = REAL( fill_value, KIND = wp )
885	ENDIF
886	to_be_resorted => km_av
887	ENDIF
888
889	CASE ( 'dummy3' ) !> @todo remove later
890	IF ( av == 0 ) THEN
891	to_be_resorted => dummy3
892	ENDIF
893
894	CASE ( 'diss1' ) !> @todo remove later
895	IF ( av == 0 ) THEN
896	to_be_resorted => dummy1
897	ENDIF
898
899	CASE ( 'diss2' ) !> @todo remove later
900	IF ( av == 0 ) THEN
901	to_be_resorted => dummy2
902	ENDIF
903
904	CASE ( 'diss_prod1' ) !> @todo remove later
905	IF ( av == 0 ) THEN
906	to_be_resorted => diss_prod1
907	ENDIF
908
909	CASE ( 'diss_adve1' ) !> @todo remove later
910	IF ( av == 0 ) THEN
911	to_be_resorted => diss_adve1
912	ENDIF
913
914	CASE ( 'diss_diff1' ) !> @todo remove later
915	IF ( av == 0 ) THEN
916	to_be_resorted => diss_diff1
917	ENDIF
918
919	CASE ( 'diss_prod2' ) !> @todo remove later
920	IF ( av == 0 ) THEN
921	to_be_resorted => diss_prod2
922	ENDIF
923
924	CASE ( 'diss_adve2' ) !> @todo remove later
925	IF ( av == 0 ) THEN
926	to_be_resorted => diss_adve2
927	ENDIF
928
929	CASE ( 'diss_diff2' ) !> @todo remove later
930	IF ( av == 0 ) THEN
931	to_be_resorted => diss_diff2
932	ENDIF
933
934	CASE ( 'diss_prod3' ) !> @todo remove later
935	IF ( av == 0 ) THEN
936	to_be_resorted => diss_prod3
937	ENDIF
938
939	CASE ( 'diss_adve3' ) !> @todo remove later
940	IF ( av == 0 ) THEN
941	to_be_resorted => diss_adve3
942	ENDIF
943
944	CASE ( 'diss_diff3' ) !> @todo remove later
945	IF ( av == 0 ) THEN
946	to_be_resorted => diss_diff3
947	ENDIF
948
949	CASE DEFAULT
950	found = .FALSE.
951
952	END SELECT
953
954
955	IF ( found .AND. .NOT. resorted ) THEN
956	DO i = nxl, nxr
957	DO j = nys, nyn
958	DO k = nzb_do, nzt_do
959	local_pf(i,j,k) = MERGE( &
960	to_be_resorted(k,j,i), &
961	REAL( fill_value, KIND = wp ), &
962	BTEST( wall_flags_0(k,j,i), flag_nr ) )
963	ENDDO
964	ENDDO
965	ENDDO
966	resorted = .TRUE.
967	ENDIF
968
969	END SUBROUTINE tcm_data_output_3d
970
971
972	!------------------------------------------------------------------------------!
973	! Description:
974	! ------------
975	!> Allocate arrays and assign pointers.
976	!------------------------------------------------------------------------------!
977	SUBROUTINE tcm_init_arrays
978
979	USE bulk_cloud_model_mod, &
980	ONLY: collision_turbulence
981
982	USE particle_attributes, &
983	ONLY: use_sgs_for_particles, wang_kernel
984
985	USE pmc_interface, &
986	ONLY: nested_run
987
988	IMPLICIT NONE
989
990	ALLOCATE( kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
991	ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
992
993	ALLOCATE( dummy1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) !> @todo remove later
994	ALLOCATE( dummy2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
995	ALLOCATE( dummy3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
996	ALLOCATE( diss_adve1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
997	ALLOCATE( diss_adve2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
998	ALLOCATE( diss_adve3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
999	ALLOCATE( diss_prod1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1000	ALLOCATE( diss_prod2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1001	ALLOCATE( diss_prod3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1002	ALLOCATE( diss_diff1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1003	ALLOCATE( diss_diff2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1004	ALLOCATE( diss_diff3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1005	dummy1 = 0.0_wp
1006	dummy2 = 0.0_wp
1007	dummy3 = 0.0_wp
1008	diss_adve1 = 0.0_wp
1009	diss_adve2 = 0.0_wp
1010	diss_adve3 = 0.0_wp
1011	diss_prod1 = 0.0_wp
1012	diss_prod2 = 0.0_wp
1013	diss_prod3 = 0.0_wp
1014	diss_diff1 = 0.0_wp
1015	diss_diff2 = 0.0_wp
1016	diss_diff3 = 0.0_wp
1017
1018	#if defined( __nopointer )
1019	ALLOCATE( e(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1020	ALLOCATE( e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1021	ALLOCATE( te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1022
1023	#else
1024	ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1025	ALLOCATE( e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1026	ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1027	#endif
1028	!
1029	!-- Allocate arrays required for dissipation.
1030	!-- Please note, if it is a nested run, arrays need to be allocated even if
1031	!-- they do not necessarily need to be transferred, which is attributed to
1032	!-- the design of the model coupler which allocates memory for each variable.
1033	IF ( rans_mode .OR. use_sgs_for_particles .OR. wang_kernel .OR. &
1034	collision_turbulence .OR. nested_run ) THEN
1035	#if defined( __nopointer )
1036	ALLOCATE( diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1037	IF ( rans_tke_e ) THEN
1038	ALLOCATE( diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1039	ALLOCATE( tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1040	ENDIF
1041	#else
1042	ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1043	IF ( rans_tke_e .OR. nested_run ) THEN
1044	ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1045	ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1046	ENDIF
1047	#endif
1048	ENDIF
1049
1050	#if ! defined( __nopointer )
1051	!
1052	!-- Initial assignment of pointers
1053	e => e_1; e_p => e_2; te_m => e_3
1054
1055	IF ( rans_mode .OR. use_sgs_for_particles .OR. &
1056	wang_kernel .OR. collision_turbulence .OR. nested_run ) THEN
1057	diss => diss_1
1058	IF ( rans_tke_e .OR. nested_run ) THEN
1059	diss_p => diss_2; tdiss_m => diss_3
1060	ENDIF
1061	ENDIF
1062	#endif
1063
1064	END SUBROUTINE tcm_init_arrays
1065
1066
1067	!------------------------------------------------------------------------------!
1068	! Description:
1069	! ------------
1070	!> Initialization of turbulence closure module.
1071	!------------------------------------------------------------------------------!
1072	SUBROUTINE tcm_init
1073
1074	USE control_parameters, &
1075	ONLY: bc_dirichlet_l, complex_terrain, topography
1076
1077	USE model_1d_mod, &
1078	ONLY: e1d, kh1d, km1d
1079
1080	USE surface_mod, &
1081	ONLY: get_topography_top_index_ji
1082
1083	IMPLICIT NONE
1084
1085	INTEGER(iwp) :: i !< loop index
1086	INTEGER(iwp) :: j !< loop index
1087	INTEGER(iwp) :: k !< loop index
1088	INTEGER(iwp) :: nz_s_shift !< lower shift index for scalars
1089	INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
1090
1091	!
1092	!-- Initialize mixing length
1093	CALL tcm_init_mixing_length
1094	dummy3 = l_wall !> @todo remove later
1095
1096	!
1097	!-- Actions for initial runs
1098	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
1099	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
1100
1101	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
1102
1103	IF ( .NOT. rans_tke_e ) THEN
1104	!
1105	!-- Transfer initial profiles to the arrays of the 3D model
1106	DO i = nxlg, nxrg
1107	DO j = nysg, nyng
1108	e(:,j,i) = e1d
1109	kh(:,j,i) = kh1d
1110	km(:,j,i) = km1d
1111	ENDDO
1112	ENDDO
1113
1114	IF ( constant_diffusion ) THEN
1115	e = 0.0_wp
1116	ENDIF
1117
1118	ELSE
1119	!
1120	!-- In case of TKE-e closure in RANS mode, do not use e, diss, and km
1121	!-- profiles from 1D model. Instead, initialize with constant profiles
1122	IF ( constant_diffusion ) THEN
1123	km = km_constant
1124	kh = km / prandtl_number
1125	e = 0.0_wp
1126	ELSEIF ( e_init > 0.0_wp ) THEN
1127	DO i = nxlg, nxrg
1128	DO j = nysg, nyng
1129	DO k = nzb+1, nzt
1130	km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1131	ENDDO
1132	ENDDO
1133	ENDDO
1134	km(nzb,:,:) = km(nzb+1,:,:)
1135	km(nzt+1,:,:) = km(nzt,:,:)
1136	kh = km / prandtl_number
1137	e = e_init
1138	ELSE
1139	IF ( .NOT. ocean_mode ) THEN
1140	kh = 0.01_wp ! there must exist an initial diffusion, because
1141	km = 0.01_wp ! otherwise no TKE would be produced by the
1142	! production terms, as long as not yet
1143	! e = (u/cm)*2 at k=nzb+1
1144	ELSE
1145	kh = 0.00001_wp
1146	km = 0.00001_wp
1147	ENDIF
1148	e = 0.0_wp
1149	ENDIF
1150
1151	DO i = nxlg, nxrg
1152	DO j = nysg, nyng
1153	DO k = nzb+1, nzt
1154	diss(k,j,i) = c_0*4 e(k,j,i)**2 / km(k,j,i)
1155	ENDDO
1156	ENDDO
1157	ENDDO
1158	diss(nzb,:,:) = diss(nzb+1,:,:)
1159	diss(nzt+1,:,:) = diss(nzt,:,:)
1160
1161	ENDIF
1162
1163	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 .OR. &
1164	INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN
1165
1166	IF ( constant_diffusion ) THEN
1167	km = km_constant
1168	kh = km / prandtl_number
1169	e = 0.0_wp
1170	ELSEIF ( e_init > 0.0_wp ) THEN
1171	DO i = nxlg, nxrg
1172	DO j = nysg, nyng
1173	DO k = nzb+1, nzt
1174	km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1175	ENDDO
1176	ENDDO
1177	ENDDO
1178	km(nzb,:,:) = km(nzb+1,:,:)
1179	km(nzt+1,:,:) = km(nzt,:,:)
1180	kh = km / prandtl_number
1181	e = e_init
1182	ELSE
1183	IF ( .NOT. ocean_mode ) THEN
1184	kh = 0.01_wp ! there must exist an initial diffusion, because
1185	km = 0.01_wp ! otherwise no TKE would be produced by the
1186	! production terms, as long as not yet
1187	! e = (u/cm)*2 at k=nzb+1
1188	ELSE
1189	kh = 0.00001_wp
1190	km = 0.00001_wp
1191	ENDIF
1192	e = 0.0_wp
1193	ENDIF
1194
1195	IF ( rans_tke_e ) THEN
1196	DO i = nxlg, nxrg
1197	DO j = nysg, nyng
1198	DO k = nzb+1, nzt
1199	diss(k,j,i) = c_0*4 e(k,j,i)**2 / km(k,j,i)
1200	ENDDO
1201	ENDDO
1202	ENDDO
1203	diss(nzb,:,:) = diss(nzb+1,:,:)
1204	diss(nzt+1,:,:) = diss(nzt,:,:)
1205	ENDIF
1206
1207	ENDIF
1208	!
1209	!-- Store initial profiles for output purposes etc.
1210	hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
1211	hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
1212	!
1213	!-- Initialize old and new time levels.
1214	te_m = 0.0_wp
1215	e_p = e
1216	IF ( rans_tke_e ) THEN
1217	tdiss_m = 0.0_wp
1218	diss_p = diss
1219	ENDIF
1220
1221	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &
1222	TRIM( initializing_actions ) == 'cyclic_fill' ) &
1223	THEN
1224
1225	!
1226	!-- In case of complex terrain and cyclic fill method as initialization,
1227	!-- shift initial data in the vertical direction for each point in the
1228	!-- x-y-plane depending on local surface height
1229	IF ( complex_terrain .AND. &
1230	TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
1231	DO i = nxlg, nxrg
1232	DO j = nysg, nyng
1233	nz_s_shift = get_topography_top_index_ji( j, i, 's' )
1234
1235	e(nz_s_shift:nzt+1,j,i) = e(0:nzt+1-nz_s_shift,j,i)
1236	km(nz_s_shift:nzt+1,j,i) = km(0:nzt+1-nz_s_shift,j,i)
1237	kh(nz_s_shift:nzt+1,j,i) = kh(0:nzt+1-nz_s_shift,j,i)
1238	ENDDO
1239	ENDDO
1240	IF ( rans_tke_e ) THEN
1241	DO i = nxlg, nxrg
1242	DO j = nysg, nyng
1243	nz_s_shift = get_topography_top_index_ji( j, i, 's' )
1244
1245	diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
1246	ENDDO
1247	ENDDO
1248	ENDIF
1249	ENDIF
1250
1251	!
1252	!-- Initialization of the turbulence recycling method
1253	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
1254	turbulent_inflow ) THEN
1255	mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e
1256	!
1257	!-- In case of complex terrain, determine vertical displacement at inflow
1258	!-- boundary and adjust mean inflow profiles
1259	IF ( complex_terrain ) THEN
1260	IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. &
1261	nysg <= 0 .AND. nyng >= 0 ) THEN
1262	nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
1263	ELSE
1264	nz_s_shift_l = 0
1265	ENDIF
1266	#if defined( __parallel )
1267	CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, &
1268	MPI_MAX, comm2d, ierr)
1269	#else
1270	nz_s_shift = nz_s_shift_l
1271	#endif
1272	mean_inflow_profiles(nz_s_shift:nzt+1,5) = &
1273	hom_sum(0:nzt+1-nz_s_shift,8,0) ! e
1274	ENDIF
1275	!
1276	!-- Use these mean profiles at the inflow (provided that Dirichlet
1277	!-- conditions are used)
1278	IF ( bc_dirichlet_l ) THEN
1279	DO j = nysg, nyng
1280	DO k = nzb, nzt+1
1281	e(k,j,nxlg:-1) = mean_inflow_profiles(k,5)
1282	ENDDO
1283	ENDDO
1284	ENDIF
1285	ENDIF
1286	!
1287	!-- Inside buildings set TKE back to zero
1288	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
1289	topography /= 'flat' ) THEN
1290	!
1291	!-- Inside buildings set TKE back to zero.
1292	!-- Other scalars (km, kh,...) are ignored at present,
1293	!-- maybe revise later.
1294	DO i = nxlg, nxrg
1295	DO j = nysg, nyng
1296	DO k = nzb, nzt
1297	e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, &
1298	BTEST( wall_flags_0(k,j,i), 0 ) )
1299	ENDDO
1300	ENDDO
1301	ENDDO
1302
1303	IF ( rans_tke_e ) THEN
1304	DO i = nxlg, nxrg
1305	DO j = nysg, nyng
1306	DO k = nzb, nzt
1307	diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp, &
1308	BTEST( wall_flags_0(k,j,i), 0 ) )
1309	ENDDO
1310	ENDDO
1311	ENDDO
1312	ENDIF
1313	ENDIF
1314	!
1315	!-- Initialize new time levels (only done in order to set boundary values
1316	!-- including ghost points)
1317	e_p = e
1318	!
1319	!-- Allthough tendency arrays are set in prognostic_equations, they have
1320	!-- to be predefined here because there they are used (but multiplied with 0)
1321	!-- before they are set.
1322	te_m = 0.0_wp
1323
1324	IF ( rans_tke_e ) THEN
1325	diss_p = diss
1326	tdiss_m = 0.0_wp
1327	ENDIF
1328
1329	ENDIF
1330
1331	END SUBROUTINE tcm_init
1332
1333
1334	! Description:
1335	! -----------------------------------------------------------------------------!
1336	!> Pre-computation of grid-dependent and near-wall mixing length.
1337	!> @todo consider walls in horizontal direction at a distance further than a
1338	!> single grid point (RANS mode)
1339	!------------------------------------------------------------------------------!
1340	SUBROUTINE tcm_init_mixing_length
1341
1342	USE arrays_3d, &
1343	ONLY: dzw, ug, vg, zu, zw
1344
1345	USE control_parameters, &
1346	ONLY: bc_lr_cyc, bc_ns_cyc, f, message_string, wall_adjustment_factor
1347
1348	USE grid_variables, &
1349	ONLY: dx, dy
1350
1351	USE indices, &
1352	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, &
1353	nzt, wall_flags_0
1354
1355	USE kinds
1356
1357
1358	IMPLICIT NONE
1359
1360	INTEGER(iwp) :: dist_dx !< found distance devided by dx
1361	INTEGER(iwp) :: i !< index variable along x
1362	INTEGER(iwp) :: ii !< index variable along x
1363	INTEGER(iwp) :: j !< index variable along y
1364	INTEGER(iwp) :: k !< index variable along z
1365	INTEGER(iwp) :: k_max_topo = 0 !< index of maximum topography height
1366	INTEGER(iwp) :: kk !< index variable along z
1367	INTEGER(iwp) :: rad_i !< search radius in grid points along x
1368	INTEGER(iwp) :: rad_i_l !< possible search radius to the left
1369	INTEGER(iwp) :: rad_i_r !< possible search radius to the right
1370	INTEGER(iwp) :: rad_j !< search radius in grid points along y
1371	INTEGER(iwp) :: rad_j_n !< possible search radius to north
1372	INTEGER(iwp) :: rad_j_s !< possible search radius to south
1373	INTEGER(iwp) :: rad_k !< search radius in grid points along z
1374	INTEGER(iwp) :: rad_k_b !< search radius in grid points along negative z
1375	INTEGER(iwp) :: rad_k_t !< search radius in grid points along positive z
1376
1377	INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
1378
1379	INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
1380
1381	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_dummy !< dummy array required for MPI_ALLREDUCE command
1382
1383	REAL(wp) :: radius !< search radius in meter
1384
1385	ALLOCATE( l_grid(1:nzt) )
1386	ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1387	!
1388	!-- Initialize the mixing length in case of an LES-simulation
1389	IF ( .NOT. rans_mode ) THEN
1390	!
1391	!-- Compute the grid-dependent mixing length.
1392	DO k = 1, nzt
1393	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp
1394	ENDDO
1395	!
1396	!-- Initialize near-wall mixing length l_wall only in the vertical direction
1397	!-- for the moment, multiplication with wall_adjustment_factor further below
1398	l_wall(nzb,:,:) = l_grid(1)
1399	DO k = nzb+1, nzt
1400	l_wall(k,:,:) = l_grid(k)
1401	ENDDO
1402	l_wall(nzt+1,:,:) = l_grid(nzt)
1403
1404	DO k = 1, nzt
1405	IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. &
1406	l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN
1407	WRITE( message_string, * ) 'grid anisotropy exceeds ', &
1408	'threshold given by only local', &
1409	' &horizontal reduction of near_wall ', &
1410	'mixing length l_wall', &
1411	' &starting from height level k = ', k, &
1412	'.'
1413	CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
1414	EXIT
1415	ENDIF
1416	ENDDO
1417	!
1418	!-- In case of topography: limit near-wall mixing length l_wall further:
1419	!-- Go through all points of the subdomain one by one and look for the closest
1420	!-- surface.
1421	!-- Is this correct in the ocean case?
1422	DO i = nxl, nxr
1423	DO j = nys, nyn
1424	DO k = nzb+1, nzt
1425	!
1426	!-- Check if current gridpoint belongs to the atmosphere
1427	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1428	!
1429	!-- Check for neighbouring grid-points.
1430	!-- Vertical distance, down
1431	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) &
1432	l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
1433	!
1434	!-- Vertical distance, up
1435	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
1436	l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
1437	!
1438	!-- y-distance
1439	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) .OR. &
1440	.NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) &
1441	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
1442	!
1443	!-- x-distance
1444	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) .OR. &
1445	.NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) &
1446	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
1447	!
1448	!-- yz-distance (vertical edges, down)
1449	IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 ) .OR. &
1450	.NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 ) ) &
1451	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1452	SQRT( 0.25_wp * dy**2 + &
1453	( zu(k) - zw(k-1) )**2 ) )
1454	!
1455	!-- yz-distance (vertical edges, up)
1456	IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 ) .OR. &
1457	.NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 ) ) &
1458	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1459	SQRT( 0.25_wp * dy**2 + &
1460	( zw(k) - zu(k) )**2 ) )
1461	!
1462	!-- xz-distance (vertical edges, down)
1463	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 ) .OR. &
1464	.NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 ) ) &
1465	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1466	SQRT( 0.25_wp * dx**2 + &
1467	( zu(k) - zw(k-1) )**2 ) )
1468	!
1469	!-- xz-distance (vertical edges, up)
1470	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 ) .OR. &
1471	.NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 ) ) &
1472	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1473	SQRT( 0.25_wp * dx**2 + &
1474	( zw(k) - zu(k) )**2 ) )
1475	!
1476	!-- xy-distance (horizontal edges)
1477	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 ) .OR. &
1478	.NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 ) .OR. &
1479	.NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 ) .OR. &
1480	.NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) ) &
1481	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1482	SQRT( 0.25_wp * ( dx2 + dy2 ) ) )
1483	!
1484	!-- xyz distance (vertical and horizontal edges, down)
1485	IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 ) .OR. &
1486	.NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 ) .OR. &
1487	.NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 ) .OR. &
1488	.NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) ) &
1489	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1490	SQRT( 0.25_wp * ( dx2 + dy2 ) &
1491	+ ( zu(k) - zw(k-1) )**2 ) )
1492	!
1493	!-- xyz distance (vertical and horizontal edges, up)
1494	IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 ) .OR. &
1495	.NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 ) .OR. &
1496	.NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 ) .OR. &
1497	.NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) ) &
1498	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1499	SQRT( 0.25_wp * ( dx2 + dy2 ) &
1500	+ ( zw(k) - zu(k) )**2 ) )
1501
1502	ENDIF
1503	ENDDO
1504	ENDDO
1505	ENDDO
1506
1507	ELSE
1508	!
1509	!-- Initialize the mixing length in case of a RANS simulation
1510	ALLOCATE( l_black(nzb:nzt+1) )
1511
1512	!
1513	!-- Calculate mixing length according to Blackadar (1962)
1514	IF ( f /= 0.0_wp ) THEN
1515	l_max = 2.7E-4_wp * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
1516	ABS( f ) + 1.0E-10_wp
1517	ELSE
1518	l_max = 30.0_wp
1519	ENDIF
1520
1521	DO k = nzb, nzt
1522	l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max )
1523	ENDDO
1524
1525	l_black(nzt+1) = l_black(nzt)
1526
1527	!
1528	!-- Get height level of highest topography within local subdomain
1529	DO i = nxlg, nxrg
1530	DO j = nysg, nyng
1531	DO k = nzb+1, nzt-1
1532	IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND. &
1533	k > k_max_topo ) &
1534	k_max_topo = k
1535	ENDDO
1536	ENDDO
1537	ENDDO
1538
1539	l_wall(nzb,:,:) = l_black(nzb)
1540	l_wall(nzt+1,:,:) = l_black(nzt+1)
1541	!
1542	!-- Limit mixing length to either nearest wall or Blackadar mixing length.
1543	!-- For that, analyze each grid point (i/j/k) ("analysed grid point") and
1544	!-- search within its vicinity for the shortest distance to a wall by cal-
1545	!-- culating the distance between the analysed grid point and the "viewed
1546	!-- grid point" if it contains a wall (belongs to topography).
1547	DO k = nzb+1, nzt
1548
1549	radius = l_black(k) ! radius within walls are searched
1550	!
1551	!-- Set l_wall to its default maximum value (l_back)
1552	l_wall(k,:,:) = radius
1553
1554	!
1555	!-- Compute search radius as number of grid points in all directions
1556	rad_i = CEILING( radius / dx )
1557	rad_j = CEILING( radius / dy )
1558
1559	DO kk = 0, nzt-k
1560	rad_k_t = kk
1561	!
1562	!-- Limit upward search radius to height of maximum topography
1563	IF ( zu(k+kk)-zu(k) >= radius .OR. k+kk >= k_max_topo ) EXIT
1564	ENDDO
1565
1566	DO kk = 0, k
1567	rad_k_b = kk
1568	IF ( zu(k)-zu(k-kk) >= radius ) EXIT
1569	ENDDO
1570
1571	!
1572	!-- Get maximum vertical radius; necessary for defining arrays
1573	rad_k = MAX( rad_k_b, rad_k_t )
1574	!
1575	!-- When analysed grid point lies above maximum topography, set search
1576	!-- radius to 0 if the distance between the analysed grid point and max
1577	!-- topography height is larger than the maximum search radius
1578	IF ( zu(k-rad_k_b) > zu(k_max_topo) ) rad_k_b = 0
1579	!
1580	!-- Search within vicinity only if the vertical search radius is >0
1581	IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 ) THEN
1582
1583	!> @note shape of vicinity is larger in z direction
1584	!> Shape of vicinity is two grid points larger than actual search
1585	!> radius in vertical direction. The first and last grid point is
1586	!> always set to 1 to asure correct detection of topography. See
1587	!> function "shortest_distance" for details.
1588	!> 2018-03-16, gronemeier
1589	ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
1590	ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
1591
1592	vicinity = 1
1593
1594	DO i = nxl, nxr
1595	DO j = nys, nyn
1596	!
1597	!-- Start search only if (i/j/k) belongs to atmosphere
1598	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1599	!
1600	!-- Reset topography within vicinity
1601	vicinity(-rad_k:rad_k,:,:) = 0
1602	!
1603	!-- Copy area surrounding analysed grid point into vicinity.
1604	!-- First, limit size of data copied to vicinity by the domain
1605	!-- border
1606	!> @note limit copied area to 1 grid point in hor. dir.
1607	!> Ignore walls in horizontal direction which are
1608	!> further away than a single grid point. This allows to
1609	!> only search within local subdomain without the need
1610	!> of global topography information.
1611	!> The error made by this assumption are acceptable at
1612	!> the moment.
1613	!> 2018-10-01, gronemeier
1614	rad_i_l = MIN( 1, rad_i, i )
1615	rad_i_r = MIN( 1, rad_i, nx-i )
1616
1617	rad_j_s = MIN( 1, rad_j, j )
1618	rad_j_n = MIN( 1, rad_j, ny-j )
1619
1620	CALL copy_into_vicinity( k, j, i, &
1621	-rad_k_b, rad_k_t, &
1622	-rad_j_s, rad_j_n, &
1623	-rad_i_l, rad_i_r )
1624	!> @note in case of cyclic boundaries, those parts of the
1625	!> topography which lies beyond the domain borders but
1626	!> still within the search radius still needs to be
1627	!> copied into 'vicinity'. As the effective search
1628	!> radius is limited to 1 at the moment, no further
1629	!> copying is needed. Old implementation (prior to
1630	!> 2018-10-01) had this covered but used a global array.
1631	!> 2018-10-01, gronemeier
1632
1633	!
1634	!-- Search for walls only if there is any within vicinity
1635	IF ( MAXVAL( vicinity(-rad_k:rad_k,:,:) ) /= 0 ) THEN
1636	!
1637	!-- Search within first half (positive x)
1638	dist_dx = rad_i
1639	DO ii = 0, dist_dx
1640	!
1641	!-- Search along vertical direction only if below
1642	!-- maximum topography
1643	IF ( rad_k_t > 0 ) THEN
1644	!
1645	!-- Search for walls within octant (+++)
1646	vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1647	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1648	shortest_distance( vic_yz, .TRUE., ii ) )
1649	!
1650	!-- Search for walls within octant (+-+)
1651	!-- Switch order of array so that the analysed grid
1652	!-- point is always located at (0/0) (required by
1653	!-- shortest_distance").
1654	vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1655	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1656	shortest_distance( vic_yz, .TRUE., ii ) )
1657
1658	ENDIF
1659	!
1660	!-- Search for walls within octant (+--)
1661	vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1662	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1663	shortest_distance( vic_yz, .FALSE., ii ) )
1664	!
1665	!-- Search for walls within octant (++-)
1666	vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1667	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1668	shortest_distance( vic_yz, .FALSE., ii ) )
1669	!
1670	!-- Reduce search along x by already found distance
1671	dist_dx = CEILING( l_wall(k,j,i) / dx )
1672
1673	ENDDO
1674	!
1675	!- Search within second half (negative x)
1676	DO ii = 0, -dist_dx, -1
1677	!
1678	!-- Search along vertical direction only if below
1679	!-- maximum topography
1680	IF ( rad_k_t > 0 ) THEN
1681	!
1682	!-- Search for walls within octant (-++)
1683	vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1684	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1685	shortest_distance( vic_yz, .TRUE., -ii ) )
1686	!
1687	!-- Search for walls within octant (--+)
1688	!-- Switch order of array so that the analysed grid
1689	!-- point is always located at (0/0) (required by
1690	!-- shortest_distance").
1691	vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1692	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1693	shortest_distance( vic_yz, .TRUE., -ii ) )
1694
1695	ENDIF
1696	!
1697	!-- Search for walls within octant (---)
1698	vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1699	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1700	shortest_distance( vic_yz, .FALSE., -ii ) )
1701	!
1702	!-- Search for walls within octant (-+-)
1703	vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1704	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1705	shortest_distance( vic_yz, .FALSE., -ii ) )
1706	!
1707	!-- Reduce search along x by already found distance
1708	dist_dx = CEILING( l_wall(k,j,i) / dx )
1709
1710	ENDDO
1711
1712	ENDIF !Check for any walls within vicinity
1713
1714	ELSE !Check if (i,j,k) belongs to atmosphere
1715
1716	l_wall(k,j,i) = l_black(k)
1717
1718	ENDIF
1719
1720	ENDDO !j loop
1721	ENDDO !i loop
1722
1723	DEALLOCATE( vicinity )
1724	DEALLOCATE( vic_yz )
1725
1726	ENDIF !check vertical size of vicinity
1727
1728	ENDDO !k loop
1729
1730	ENDIF !LES or RANS mode
1731
1732	!
1733	!-- Set lateral boundary conditions for l_wall
1734	CALL exchange_horiz( l_wall, nbgp )
1735
1736	CONTAINS
1737	!------------------------------------------------------------------------------!
1738	! Description:
1739	! ------------
1740	!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and
1741	!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
1742	!> closest to the origin (0/0) of 'array'.
1743	!------------------------------------------------------------------------------!
1744	REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
1745
1746	IMPLICIT NONE
1747
1748	LOGICAL, INTENT(IN) :: orientation !< flag if array represents an array oriented upwards (true) or downwards (false)
1749
1750	INTEGER(iwp), INTENT(IN) :: pos_i !< x position of the yz-plane 'array'
1751
1752	INTEGER(iwp) :: a !< loop index
1753	INTEGER(iwp) :: b !< loop index
1754	INTEGER(iwp) :: jj !< loop index
1755
1756	INTEGER(KIND=1) :: maximum !< maximum of array along z dimension
1757
1758	INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
1759
1760	INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i
1761
1762	!
1763	!-- Get coordinate of first maximum along vertical dimension
1764	!-- at each y position of array (similar to function maxloc but more stable).
1765	DO a = 0, rad_j
1766	loc_k(a) = rad_k+1
1767	maximum = MAXVAL( array(:,a) )
1768	DO b = 0, rad_k+1
1769	IF ( array(b,a) == maximum ) THEN
1770	loc_k(a) = b
1771	EXIT
1772	ENDIF
1773	ENDDO
1774	ENDDO
1775	!
1776	!-- Set distance to the default maximum value (=search radius)
1777	shortest_distance = radius
1778	!
1779	!-- Calculate distance between position (0/0/0) and
1780	!-- position (pos_i/jj/loc(jj)) and only save the shortest distance.
1781	IF ( orientation ) THEN !if array is oriented upwards
1782	DO jj = 0, rad_j
1783	shortest_distance = &
1784	MIN( shortest_distance, &
1785	SQRT( MAX(REAL(pos_i, KIND=wp)dx-0.5_wpdx, 0.0_wp)**2 &
1786	+ MAX(REAL(jj, KIND=wp)dy-0.5_wpdy, 0.0_wp)**2 &
1787	+ MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2 &
1788	) &
1789	)
1790	ENDDO
1791	ELSE !if array is oriented downwards
1792	!> @note MAX within zw required to circumvent error at domain border
1793	!> At the domain border, if non-cyclic boundary is present, the
1794	!> index for zw could be -1, which will be errorneous (zw(-1) does
1795	!> not exist). The MAX function limits the index to be at least 0.
1796	DO jj = 0, rad_j
1797	shortest_distance = &
1798	MIN( shortest_distance, &
1799	SQRT( MAX(REAL(pos_i, KIND=wp)dx-0.5_wpdx, 0.0_wp)**2 &
1800	+ MAX(REAL(jj, KIND=wp)dy-0.5_wpdy, 0.0_wp)**2 &
1801	+ MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2 &
1802	) &
1803	)
1804	ENDDO
1805	ENDIF
1806
1807	END FUNCTION
1808
1809	!------------------------------------------------------------------------------!
1810	! Description:
1811	! ------------
1812	!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point
1813	!> (kp,jp,ip) containing the first bit of wall_flags_0 into the array
1814	!> 'vicinity'. Only copy first bit as this indicates the presence of topography.
1815	!------------------------------------------------------------------------------!
1816	SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir )
1817
1818	IMPLICIT NONE
1819
1820	INTEGER(iwp), INTENT(IN) :: il !< left loop boundary
1821	INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction
1822	INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary
1823	INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary
1824	INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction
1825	INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary
1826	INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary
1827	INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction
1828	INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary
1829
1830	INTEGER(iwp) :: i !< loop index
1831	INTEGER(iwp) :: j !< loop index
1832	INTEGER(iwp) :: k !< loop index
1833
1834	DO i = il, ir
1835	DO j = js, jn
1836	DO k = kb, kt
1837	vicinity(k,j,i) = MERGE( 0, 1, &
1838	BTEST( wall_flags_0(kp+k,jp+j,ip+i), 0 ) )
1839	ENDDO
1840	ENDDO
1841	ENDDO
1842
1843	END SUBROUTINE copy_into_vicinity
1844
1845	END SUBROUTINE tcm_init_mixing_length
1846
1847
1848	!------------------------------------------------------------------------------!
1849	! Description:
1850	! ------------
1851	!> Initialize virtual velocities used later in production_e.
1852	!------------------------------------------------------------------------------!
1853	SUBROUTINE production_e_init
1854
1855	USE arrays_3d, &
1856	ONLY: drho_air_zw, zu
1857
1858	USE control_parameters, &
1859	ONLY: constant_flux_layer
1860
1861	USE surface_layer_fluxes_mod, &
1862	ONLY: phi_m
1863
1864	USE surface_mod, &
1865	ONLY : surf_def_h, surf_lsm_h, surf_usm_h
1866
1867	IMPLICIT NONE
1868
1869	INTEGER(iwp) :: i !< grid index x-direction
1870	INTEGER(iwp) :: j !< grid index y-direction
1871	INTEGER(iwp) :: k !< grid index z-direction
1872	INTEGER(iwp) :: m !< running index surface elements
1873
1874	REAL(wp) :: km_sfc !< diffusion coefficient, used to compute virtual velocities
1875
1876	IF ( constant_flux_layer ) THEN
1877	!
1878	!-- Calculate a virtual velocity at the surface in a way that the
1879	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1880	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1881	!-- production term at k=1 (see production_e_ij).
1882	!-- The velocity gradient has to be limited in case of too small km
1883	!-- (otherwise the timestep may be significantly reduced by large
1884	!-- surface winds).
1885	!-- not available in case of non-cyclic boundary conditions.
1886	!-- Default surfaces, upward-facing
1887	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1888	DO m = 1, surf_def_h(0)%ns
1889
1890	i = surf_def_h(0)%i(m)
1891	j = surf_def_h(0)%j(m)
1892	k = surf_def_h(0)%k(m)
1893	!
1894	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1895	!-- and km are not on the same grid. Actually, a further
1896	!-- interpolation of km onto the u/v-grid is necessary. However, the
1897	!-- effect of this error is negligible.
1898	km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / &
1899	phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
1900
1901	surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * &
1902	drho_air_zw(k-1) * &
1903	( zu(k+1) - zu(k-1) ) / &
1904	( km_sfc + 1.0E-20_wp )
1905	surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * &
1906	drho_air_zw(k-1) * &
1907	( zu(k+1) - zu(k-1) ) / &
1908	( km_sfc + 1.0E-20_wp )
1909
1910	IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > &
1911	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1912	) surf_def_h(0)%u_0(m) = u(k-1,j,i)
1913
1914	IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > &
1915	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1916	) surf_def_h(0)%v_0(m) = v(k-1,j,i)
1917
1918	ENDDO
1919	!
1920	!-- Default surfaces, downward-facing surfaces
1921	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1922	DO m = 1, surf_def_h(1)%ns
1923
1924	i = surf_def_h(1)%i(m)
1925	j = surf_def_h(1)%j(m)
1926	k = surf_def_h(1)%k(m)
1927	!
1928	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1929	!-- and km are not on the same grid. Actually, a further
1930	!-- interpolation of km onto the u/v-grid is necessary. However, the
1931	!-- effect of this error is negligible.
1932	surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * &
1933	drho_air_zw(k-1) * &
1934	( zu(k+1) - zu(k-1) ) / &
1935	( km(k,j,i) + 1.0E-20_wp )
1936	surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * &
1937	drho_air_zw(k-1) * &
1938	( zu(k+1) - zu(k-1) ) / &
1939	( km(k,j,i) + 1.0E-20_wp )
1940
1941	IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > &
1942	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1943	) surf_def_h(1)%u_0(m) = u(k+1,j,i)
1944
1945	IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > &
1946	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1947	) surf_def_h(1)%v_0(m) = v(k+1,j,i)
1948
1949	ENDDO
1950	!
1951	!-- Natural surfaces, upward-facing
1952	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1953	DO m = 1, surf_lsm_h%ns
1954
1955	i = surf_lsm_h%i(m)
1956	j = surf_lsm_h%j(m)
1957	k = surf_lsm_h%k(m)
1958	!
1959	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1960	!-- and km are not on the same grid. Actually, a further
1961	!-- interpolation of km onto the u/v-grid is necessary. However, the
1962	!-- effect of this error is negligible.
1963	km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / &
1964	phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
1965
1966	surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * &
1967	drho_air_zw(k-1) * &
1968	( zu(k+1) - zu(k-1) ) / &
1969	( km_sfc + 1.0E-20_wp )
1970	surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * &
1971	drho_air_zw(k-1) * &
1972	( zu(k+1) - zu(k-1) ) / &
1973	( km_sfc + 1.0E-20_wp )
1974
1975	IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > &
1976	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1977	) surf_lsm_h%u_0(m) = u(k-1,j,i)
1978
1979	IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > &
1980	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1981	) surf_lsm_h%v_0(m) = v(k-1,j,i)
1982
1983	ENDDO
1984	!
1985	!-- Urban surfaces, upward-facing
1986	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1987	DO m = 1, surf_usm_h%ns
1988
1989	i = surf_usm_h%i(m)
1990	j = surf_usm_h%j(m)
1991	k = surf_usm_h%k(m)
1992	!
1993	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1994	!-- and km are not on the same grid. Actually, a further
1995	!-- interpolation of km onto the u/v-grid is necessary. However, the
1996	!-- effect of this error is negligible.
1997	km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / &
1998	phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
1999
2000	surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * &
2001	drho_air_zw(k-1) * &
2002	( zu(k+1) - zu(k-1) ) / &
2003	( km_sfc + 1.0E-20_wp )
2004	surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * &
2005	drho_air_zw(k-1) * &
2006	( zu(k+1) - zu(k-1) ) / &
2007	( km_sfc + 1.0E-20_wp )
2008
2009	IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > &
2010	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
2011	) surf_usm_h%u_0(m) = u(k-1,j,i)
2012
2013	IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > &
2014	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
2015	) surf_usm_h%v_0(m) = v(k-1,j,i)
2016
2017	ENDDO
2018
2019	ENDIF
2020
2021	END SUBROUTINE production_e_init
2022
2023
2024	!------------------------------------------------------------------------------!
2025	! Description:
2026	! ------------
2027	!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2028	!> Vector-optimized version
2029	!------------------------------------------------------------------------------!
2030	SUBROUTINE tcm_prognostic_equations
2031
2032	USE arrays_3d, &
2033	ONLY: ddzu
2034
2035	USE control_parameters, &
2036	ONLY: f, scalar_advec, tsc
2037
2038	USE surface_mod, &
2039	ONLY : surf_def_h
2040
2041	IMPLICIT NONE
2042
2043	INTEGER(iwp) :: i !< loop index
2044	INTEGER(iwp) :: j !< loop index
2045	INTEGER(iwp) :: k !< loop index
2046	INTEGER(iwp) :: m !< loop index
2047	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
2048	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
2049
2050	REAL(wp) :: sbt !< wheighting factor for sub-time step
2051
2052	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: advec !< advection term of TKE tendency
2053	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: produc !< production term of TKE tendency
2054
2055	!
2056	!-- If required, compute prognostic equation for turbulent kinetic
2057	!-- energy (TKE)
2058	IF ( .NOT. constant_diffusion ) THEN
2059
2060	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
2061
2062	sbt = tsc(2)
2063	IF ( .NOT. use_upstream_for_tke ) THEN
2064	IF ( scalar_advec == 'bc-scheme' ) THEN
2065
2066	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2067	!
2068	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2069	sbt = 1.0_wp
2070	ENDIF
2071	tend = 0.0_wp
2072	CALL advec_s_bc( e, 'e' )
2073
2074	ENDIF
2075	ENDIF
2076
2077	!
2078	!-- TKE-tendency terms with no communication
2079	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
2080	IF ( use_upstream_for_tke ) THEN
2081	tend = 0.0_wp
2082	CALL advec_s_up( e )
2083	ELSE
2084	tend = 0.0_wp
2085	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2086	IF ( ws_scheme_sca ) THEN
2087	CALL advec_s_ws( e, 'e' )
2088	ELSE
2089	CALL advec_s_pw( e )
2090	ENDIF
2091	ELSE
2092	CALL advec_s_up( e )
2093	ENDIF
2094	ENDIF
2095	ENDIF
2096
2097	IF ( rans_tke_e ) advec = tend
2098
2099	CALL production_e( .FALSE. )
2100
2101	!
2102	!-- Save production term for prognostic equation of TKE dissipation rate
2103	IF ( rans_tke_e ) produc = tend - advec
2104
2105	IF ( .NOT. humidity ) THEN
2106	IF ( ocean_mode ) THEN
2107	CALL diffusion_e( prho, prho_reference )
2108	ELSE
2109	CALL diffusion_e( pt, pt_reference )
2110	ENDIF
2111	ELSE
2112	CALL diffusion_e( vpt, pt_reference )
2113	ENDIF
2114
2115	!
2116	!-- Additional sink term for flows through plant canopies
2117	IF ( plant_canopy ) CALL pcm_tendency( 6 )
2118
2119	CALL user_actions( 'e-tendency' )
2120
2121	!
2122	!-- Prognostic equation for TKE.
2123	!-- Eliminate negative TKE values, which can occur due to numerical
2124	!-- reasons in the course of the integration. In such cases the old TKE
2125	!-- value is reduced by 90%.
2126	DO i = nxl, nxr
2127	DO j = nys, nyn
2128	DO k = nzb+1, nzt
2129	e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2130	tsc(3) * te_m(k,j,i) ) &
2131	) &
2132	* MERGE( 1.0_wp, 0.0_wp, &
2133	BTEST( wall_flags_0(k,j,i), 0 ) &
2134	)
2135	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
2136	ENDDO
2137	ENDDO
2138	ENDDO
2139
2140	!
2141	!-- Use special boundary condition in case of TKE-e closure
2142	!> @todo do the same for usm and lsm surfaces
2143	!> 2018-06-05, gronemeier
2144	IF ( rans_tke_e ) THEN
2145	DO i = nxl, nxr
2146	DO j = nys, nyn
2147	surf_s = surf_def_h(0)%start_index(j,i)
2148	surf_e = surf_def_h(0)%end_index(j,i)
2149	DO m = surf_s, surf_e
2150	k = surf_def_h(0)%k(m)
2151	e_p(k,j,i) = surf_def_h(0)%us(m)2 / c_02
2152	ENDDO
2153	ENDDO
2154	ENDDO
2155	ENDIF
2156
2157	!
2158	!-- Calculate tendencies for the next Runge-Kutta step
2159	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2160	IF ( intermediate_timestep_count == 1 ) THEN
2161	DO i = nxl, nxr
2162	DO j = nys, nyn
2163	DO k = nzb+1, nzt
2164	te_m(k,j,i) = tend(k,j,i)
2165	ENDDO
2166	ENDDO
2167	ENDDO
2168	ELSEIF ( intermediate_timestep_count < &
2169	intermediate_timestep_count_max ) THEN
2170	DO i = nxl, nxr
2171	DO j = nys, nyn
2172	DO k = nzb+1, nzt
2173	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2174	+ 5.3125_wp * te_m(k,j,i)
2175	ENDDO
2176	ENDDO
2177	ENDDO
2178	ENDIF
2179	ENDIF
2180
2181	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
2182
2183	ENDIF ! TKE equation
2184
2185	!
2186	!-- If required, compute prognostic equation for TKE dissipation rate
2187	IF ( rans_tke_e ) THEN
2188
2189	CALL cpu_log( log_point(33), 'diss-equation', 'start' )
2190
2191	sbt = tsc(2)
2192	IF ( .NOT. use_upstream_for_tke ) THEN
2193	IF ( scalar_advec == 'bc-scheme' ) THEN
2194
2195	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2196	!
2197	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2198	sbt = 1.0_wp
2199	ENDIF
2200	tend = 0.0_wp
2201	CALL advec_s_bc( diss, 'diss' )
2202
2203	ENDIF
2204	ENDIF
2205
2206	!
2207	!-- dissipation-tendency terms with no communication
2208	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
2209	IF ( use_upstream_for_tke ) THEN
2210	tend = 0.0_wp
2211	CALL advec_s_up( diss )
2212	ELSE
2213	tend = 0.0_wp
2214	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2215	IF ( ws_scheme_sca ) THEN
2216	CALL advec_s_ws( diss, 'diss' )
2217	ELSE
2218	CALL advec_s_pw( diss )
2219	ENDIF
2220	ELSE
2221	CALL advec_s_up( diss )
2222	ENDIF
2223	ENDIF
2224	ENDIF
2225
2226	!
2227	!-- Production of TKE dissipation rate
2228	DO i = nxl, nxr
2229	DO j = nys, nyn
2230	DO k = nzb+1, nzt
2231	! tend(k,j,i) = tend(k,j,i) + c_1 * diss(k,j,i) / ( e(k,j,i) + 1.0E-20_wp ) * produc(k)
2232	tend(k,j,i) = tend(k,j,i) + c_1 * c_0*4 f / c_4 & !> @todo needs revision
2233	/ surf_def_h(0)%us(surf_def_h(0)%start_index(j,i)) &
2234	* SQRT(e(k,j,i)) * produc(k,j,i)
2235	ENDDO
2236	ENDDO
2237	ENDDO
2238
2239	CALL diffusion_diss
2240
2241	!
2242	!-- Additional sink term for flows through plant canopies
2243	! IF ( plant_canopy ) CALL pcm_tendency( ? ) !> @query what to do with this?
2244
2245	! CALL user_actions( 'diss-tendency' ) !> @todo not yet implemented
2246
2247	!
2248	!-- Prognostic equation for TKE dissipation.
2249	!-- Eliminate negative dissipation values, which can occur due to numerical
2250	!-- reasons in the course of the integration. In such cases the old
2251	!-- dissipation value is reduced by 90%.
2252	DO i = nxl, nxr
2253	DO j = nys, nyn
2254	DO k = nzb+1, nzt
2255	diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2256	tsc(3) * tdiss_m(k,j,i) ) &
2257	) &
2258	* MERGE( 1.0_wp, 0.0_wp, &
2259	BTEST( wall_flags_0(k,j,i), 0 ) &
2260	)
2261	IF ( diss_p(k,j,i) < 0.0_wp ) &
2262	diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
2263	ENDDO
2264	ENDDO
2265	ENDDO
2266
2267	!
2268	!-- Use special boundary condition in case of TKE-e closure
2269	DO i = nxl, nxr
2270	DO j = nys, nyn
2271	surf_s = surf_def_h(0)%start_index(j,i)
2272	surf_e = surf_def_h(0)%end_index(j,i)
2273	DO m = surf_s, surf_e
2274	k = surf_def_h(0)%k(m)
2275	diss_p(k,j,i) = surf_def_h(0)%us(m)*3 / kappa ddzu(k)
2276	ENDDO
2277	ENDDO
2278	ENDDO
2279
2280	!
2281	!-- Calculate tendencies for the next Runge-Kutta step
2282	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2283	IF ( intermediate_timestep_count == 1 ) THEN
2284	DO i = nxl, nxr
2285	DO j = nys, nyn
2286	DO k = nzb+1, nzt
2287	tdiss_m(k,j,i) = tend(k,j,i)
2288	ENDDO
2289	ENDDO
2290	ENDDO
2291	ELSEIF ( intermediate_timestep_count < &
2292	intermediate_timestep_count_max ) THEN
2293	DO i = nxl, nxr
2294	DO j = nys, nyn
2295	DO k = nzb+1, nzt
2296	tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2297	+ 5.3125_wp * tdiss_m(k,j,i)
2298	ENDDO
2299	ENDDO
2300	ENDDO
2301	ENDIF
2302	ENDIF
2303
2304	CALL cpu_log( log_point(33), 'diss-equation', 'stop' )
2305
2306	ENDIF
2307
2308	END SUBROUTINE tcm_prognostic_equations
2309
2310
2311	!------------------------------------------------------------------------------!
2312	! Description:
2313	! ------------
2314	!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2315	!> Cache-optimized version
2316	!------------------------------------------------------------------------------!
2317	SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn )
2318
2319	USE arrays_3d, &
2320	ONLY: diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss, &
2321	flux_l_e, flux_s_diss, flux_s_e
2322
2323	USE control_parameters, &
2324	ONLY: tsc
2325
2326	IMPLICIT NONE
2327
2328	INTEGER(iwp) :: i !< loop index x direction
2329	INTEGER(iwp) :: i_omp !< first loop index of i-loop in prognostic_equations
2330	INTEGER(iwp) :: j !< loop index y direction
2331	INTEGER(iwp) :: k !< loop index z direction
2332	INTEGER(iwp) :: tn !< task number of openmp task
2333
2334	! INTEGER(iwp) :: pis = 32 !< debug variable, print from i=pis !> @todo remove later
2335	! INTEGER(iwp) :: pie = 32 !< debug variable, print until i=pie !> @todo remove later
2336	! INTEGER(iwp) :: pjs = 26 !< debug variable, print from j=pjs !> @todo remove later
2337	! INTEGER(iwp) :: pje = 26 !< debug variable, print until j=pje !> @todo remove later
2338	! INTEGER(iwp) :: pkb = 1 !< debug variable, print from k=pkb !> @todo remove later
2339	! INTEGER(iwp) :: pkt = 7 !< debug variable, print until k=pkt !> @todo remove later
2340
2341	REAL(wp), DIMENSION(nzb:nzt+1) :: dum_adv !< debug variable !> @todo remove later
2342	REAL(wp), DIMENSION(nzb:nzt+1) :: dum_pro !< debug variable !> @todo remove later
2343	REAL(wp), DIMENSION(nzb:nzt+1) :: dum_dif !< debug variable !> @todo remove later
2344
2345	!5555 FORMAT(A,7(1X,E12.5)) !> @todo remove later
2346
2347	!
2348	!-- If required, compute prognostic equation for turbulent kinetic
2349	!-- energy (TKE)
2350	IF ( .NOT. constant_diffusion ) THEN
2351
2352	!
2353	!-- Tendency-terms for TKE
2354	tend(:,j,i) = 0.0_wp
2355	IF ( timestep_scheme(1:5) == 'runge' &
2356	.AND. .NOT. use_upstream_for_tke ) THEN
2357	IF ( ws_scheme_sca ) THEN
2358	CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
2359	flux_l_e, diss_l_e , i_omp, tn )
2360	ELSE
2361	CALL advec_s_pw( i, j, e )
2362	ENDIF
2363	ELSE
2364	CALL advec_s_up( i, j, e )
2365	ENDIF
2366
2367	dum_adv = tend(:,j,i) !> @todo remove later
2368
2369	CALL production_e( i, j, .FALSE. )
2370
2371	dum_pro = tend(:,j,i) - dum_adv !> @todo remove later
2372
2373	IF ( .NOT. humidity ) THEN
2374	IF ( ocean_mode ) THEN
2375	CALL diffusion_e( i, j, prho, prho_reference )
2376	ELSE
2377	CALL diffusion_e( i, j, pt, pt_reference )
2378	ENDIF
2379	ELSE
2380	CALL diffusion_e( i, j, vpt, pt_reference )
2381	ENDIF
2382
2383	dum_dif = tend(:,j,i) - dum_adv - dum_pro !> @todo remove later
2384
2385	!
2386	!-- Additional sink term for flows through plant canopies
2387	IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 )
2388
2389	CALL user_actions( i, j, 'e-tendency' )
2390
2391	!
2392	!-- Prognostic equation for TKE.
2393	!-- Eliminate negative TKE values, which can occur due to numerical
2394	!-- reasons in the course of the integration. In such cases the old
2395	!-- TKE value is reduced by 90%.
2396	DO k = nzb+1, nzt
2397	e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
2398	tsc(3) * te_m(k,j,i) ) &
2399	) &
2400	* MERGE( 1.0_wp, 0.0_wp, &
2401	BTEST( wall_flags_0(k,j,i), 0 ) &
2402	)
2403	IF ( e_p(k,j,i) <= 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
2404	ENDDO
2405
2406	!
2407	!-- Calculate tendencies for the next Runge-Kutta step
2408	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2409	IF ( intermediate_timestep_count == 1 ) THEN
2410	DO k = nzb+1, nzt
2411	te_m(k,j,i) = tend(k,j,i)
2412	ENDDO
2413	ELSEIF ( intermediate_timestep_count < &
2414	intermediate_timestep_count_max ) THEN
2415	DO k = nzb+1, nzt
2416	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2417	5.3125_wp * te_m(k,j,i)
2418	ENDDO
2419	ENDIF
2420	ENDIF
2421
2422	! if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then !> @todo remove later
2423	! WRITE(9, *) '------'
2424	! WRITE(9, '(A,F8.3,1X,F8.3,1X,I2)') 't, dt, int_ts:', simulated_time, dt_3d, intermediate_timestep_count
2425	! WRITE(9, *) 'i:', i
2426	! WRITE(9, *) 'j:', j
2427	! WRITE(9, *) 'k:', pkb, ' - ', pkt
2428	! WRITE(9, *) '---'
2429	! WRITE(9, *) 'e:'
2430	! WRITE(9, 5555) 'adv :', dum_adv(pkb:pkt)
2431	! WRITE(9, 5555) 'pro :', dum_pro(pkb:pkt)
2432	! WRITE(9, 5555) 'dif :', dum_dif(pkb:pkt)
2433	! WRITE(9, 5555) 'tend:', tend(pkb:pkt,j,i)
2434	! WRITE(9, 5555) 'e_p :', e_p(pkb:pkt,j,i)
2435	! WRITE(9, 5555) 'e :', e(pkb:pkt,j,i)
2436	! FLUSH(9)
2437	! endif
2438
2439	ENDIF ! TKE equation
2440
2441	!
2442	!-- If required, compute prognostic equation for TKE dissipation rate
2443	IF ( rans_tke_e ) THEN
2444
2445	!
2446	!-- Tendency-terms for dissipation
2447	tend(:,j,i) = 0.0_wp
2448	IF ( timestep_scheme(1:5) == 'runge' &
2449	.AND. .NOT. use_upstream_for_tke ) THEN
2450	IF ( ws_scheme_sca ) THEN
2451	CALL advec_s_ws( i, j, diss, 'diss', flux_s_diss, diss_s_diss, &
2452	flux_l_diss, diss_l_diss, i_omp, tn )
2453	ELSE
2454	CALL advec_s_pw( i, j, diss )
2455	ENDIF
2456	ELSE
2457	CALL advec_s_up( i, j, diss )
2458	ENDIF
2459
2460	IF ( intermediate_timestep_count == 1 ) diss_adve1(:,j,i) = tend(:,j,i) !> @todo remove later
2461	IF ( intermediate_timestep_count == 2 ) diss_adve2(:,j,i) = tend(:,j,i)
2462	IF ( intermediate_timestep_count == 3 ) diss_adve3(:,j,i) = tend(:,j,i)
2463
2464	!
2465	!-- Production of TKE dissipation rate
2466	CALL production_e( i, j, .TRUE. )
2467
2468	IF ( intermediate_timestep_count == 1 ) diss_prod1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) !> @todo remove later
2469	IF ( intermediate_timestep_count == 2 ) diss_prod2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i)
2470	IF ( intermediate_timestep_count == 3 ) diss_prod3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i)
2471
2472	dum_pro = tend(:,j,i) - dum_adv !> @todo remove later
2473
2474	!
2475	!-- Diffusion term of TKE dissipation rate
2476	CALL diffusion_diss( i, j )
2477
2478	IF ( intermediate_timestep_count == 1 ) diss_diff1(:,j,i) = tend(:,j,i) - diss_adve1(:,j,i) - diss_prod1(:,j,i) !> @todo remove later
2479	IF ( intermediate_timestep_count == 2 ) diss_diff2(:,j,i) = tend(:,j,i) - diss_adve2(:,j,i) - diss_prod2(:,j,i)
2480	IF ( intermediate_timestep_count == 3 ) diss_diff3(:,j,i) = tend(:,j,i) - diss_adve3(:,j,i) - diss_prod3(:,j,i)
2481	! IF ( intermediate_timestep_count == 3 ) dummy3(:,j,i) = km(:,j,i)
2482
2483	dum_dif = tend(:,j,i) - dum_adv - dum_pro !> @todo remove later
2484
2485	!
2486	!-- Additional sink term for flows through plant canopies
2487	! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !> @todo not yet implemented
2488
2489	! CALL user_actions( i, j, 'diss-tendency' ) !> @todo not yet implemented
2490
2491	!
2492	!-- Prognostic equation for TKE dissipation
2493	!-- Eliminate negative dissipation values, which can occur due to
2494	!-- numerical reasons in the course of the integration. In such cases
2495	!-- the old dissipation value is reduced by 90%.
2496	DO k = nzb+1, nzt
2497	diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
2498	tsc(3) * tdiss_m(k,j,i) ) &
2499	) &
2500	* MERGE( 1.0_wp, 0.0_wp, &
2501	BTEST( wall_flags_0(k,j,i), 0 )&
2502	)
2503	ENDDO
2504
2505	!
2506	!-- Calculate tendencies for the next Runge-Kutta step
2507	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2508	IF ( intermediate_timestep_count == 1 ) THEN
2509	DO k = nzb+1, nzt
2510	tdiss_m(k,j,i) = tend(k,j,i)
2511	ENDDO
2512	ELSEIF ( intermediate_timestep_count < &
2513	intermediate_timestep_count_max ) THEN
2514	DO k = nzb+1, nzt
2515	tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2516	5.3125_wp * tdiss_m(k,j,i)
2517	ENDDO
2518	ENDIF
2519	ENDIF
2520
2521	IF ( intermediate_timestep_count == 1 ) dummy1(:,j,i) = diss_p(:,j,i) !> @todo remove later
2522	IF ( intermediate_timestep_count == 2 ) dummy2(:,j,i) = diss_p(:,j,i)
2523
2524	! if ( i >= pis .and. i <= pie .and. j >= pjs .and. j <= pje ) then !> @todo remove later
2525	! WRITE(9, *) '---'
2526	! WRITE(9, *) 'diss:'
2527	! WRITE(9, 5555) 'adv :', dum_adv(pkb:pkt)
2528	! WRITE(9, 5555) 'pro :', dum_pro(pkb:pkt)
2529	! WRITE(9, 5555) 'dif :', dum_dif(pkb:pkt)
2530	! WRITE(9, 5555) 'tend :', tend(pkb:pkt,j,i)
2531	! WRITE(9, 5555) 'diss_p:', diss_p(pkb:pkt,j,i)
2532	! WRITE(9, 5555) 'diss :', diss(pkb:pkt,j,i)
2533	! WRITE(9, *) '---'
2534	! WRITE(9, 5555) 'km :', km(pkb:pkt,j,i)
2535	! flush(9)
2536	! endif
2537
2538	ENDIF ! dissipation equation
2539
2540	END SUBROUTINE tcm_prognostic_equations_ij
2541
2542
2543	!------------------------------------------------------------------------------!
2544	! Description:
2545	! ------------
2546	!> Production terms (shear + buoyancy) of the TKE.
2547	!> Vector-optimized version
2548	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
2549	!> not considered well!
2550	!> @todo Adjust production term in case of rans_tke_e simulation
2551	!------------------------------------------------------------------------------!
2552	SUBROUTINE production_e( diss_production )
2553
2554	USE arrays_3d, &
2555	ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
2556
2557	USE control_parameters, &
2558	ONLY: cloud_droplets, constant_flux_layer, neutral, &
2559	rho_reference, use_single_reference_value, use_surface_fluxes, &
2560	use_top_fluxes
2561
2562	USE grid_variables, &
2563	ONLY: ddx, dx, ddy, dy
2564
2565	USE bulk_cloud_model_mod, &
2566	ONLY: bulk_cloud_model
2567
2568	USE surface_mod, &
2569	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
2570	surf_usm_v
2571
2572	IMPLICIT NONE
2573
2574	LOGICAL :: diss_production
2575
2576	INTEGER(iwp) :: i !< running index x-direction
2577	INTEGER(iwp) :: j !< running index y-direction
2578	INTEGER(iwp) :: k !< running index z-direction
2579	INTEGER(iwp) :: l !< running index for different surface type orientation
2580	INTEGER(iwp) :: m !< running index surface elements
2581	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
2582	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
2583	INTEGER(iwp) :: flag_nr !< number of masking flag
2584
2585	REAL(wp) :: def !<
2586	REAL(wp) :: flag !< flag to mask topography
2587	REAL(wp) :: k1 !<
2588	REAL(wp) :: k2 !<
2589	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
2590	REAL(wp) :: theta !<
2591	REAL(wp) :: temp !<
2592	REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation
2593	REAL(wp) :: usvs !< momentum flux u"v"
2594	REAL(wp) :: vsus !< momentum flux v"u"
2595	REAL(wp) :: wsus !< momentum flux w"u"
2596	REAL(wp) :: wsvs !< momentum flux w"v"
2597
2598	REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction
2599	REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction
2600	REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction
2601	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction
2602	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction
2603	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction
2604	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction
2605	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction
2606	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction
2607	REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary flux-array in z-direction
2608
2609
2610
2611	!
2612	!-- Calculate TKE production by shear. Calculate gradients at all grid
2613	!-- points first, gradients at surface-bounded grid points will be
2614	!-- overwritten further below.
2615	DO i = nxl, nxr
2616	DO j = nys, nyn
2617	DO k = nzb+1, nzt
2618
2619	dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
2620	dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
2621	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
2622	dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
2623	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
2624
2625	dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
2626	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
2627	dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
2628	dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
2629	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
2630
2631	dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
2632	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
2633	dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
2634	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
2635	dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2636
2637	ENDDO
2638
2639
2640	flag_nr = 29
2641
2642
2643	IF ( constant_flux_layer ) THEN
2644	!
2645
2646	flag_nr = 0
2647
2648	!-- Position beneath wall
2649	!-- (2) - Will allways be executed.
2650	!-- 'bottom and wall: use u_0,v_0 and wall functions'
2651	!
2652	!-- Compute gradients at north- and south-facing surfaces.
2653	!-- First, for default surfaces, then for urban surfaces.
2654	!-- Note, so far no natural vertical surfaces implemented
2655	DO l = 0, 1
2656	surf_s = surf_def_v(l)%start_index(j,i)
2657	surf_e = surf_def_v(l)%end_index(j,i)
2658	DO m = surf_s, surf_e
2659	k = surf_def_v(l)%k(m)
2660	usvs = surf_def_v(l)%mom_flux_tke(0,m)
2661	wsvs = surf_def_v(l)%mom_flux_tke(1,m)
2662
2663	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
2664	* 0.5_wp * dy
2665	!
2666	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2667	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2668	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2669	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2670	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2671	ENDDO
2672	!
2673	!-- Natural surfaces
2674	surf_s = surf_lsm_v(l)%start_index(j,i)
2675	surf_e = surf_lsm_v(l)%end_index(j,i)
2676	DO m = surf_s, surf_e
2677	k = surf_lsm_v(l)%k(m)
2678	usvs = surf_lsm_v(l)%mom_flux_tke(0,m)
2679	wsvs = surf_lsm_v(l)%mom_flux_tke(1,m)
2680
2681	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
2682	* 0.5_wp * dy
2683	!
2684	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2685	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2686	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2687	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2688	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2689	ENDDO
2690	!
2691	!-- Urban surfaces
2692	surf_s = surf_usm_v(l)%start_index(j,i)
2693	surf_e = surf_usm_v(l)%end_index(j,i)
2694	DO m = surf_s, surf_e
2695	k = surf_usm_v(l)%k(m)
2696	usvs = surf_usm_v(l)%mom_flux_tke(0,m)
2697	wsvs = surf_usm_v(l)%mom_flux_tke(1,m)
2698
2699	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
2700	* 0.5_wp * dy
2701	!
2702	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2703	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2704	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2705	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2706	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2707	ENDDO
2708	ENDDO
2709	!
2710	!-- Compute gradients at east- and west-facing walls
2711	DO l = 2, 3
2712	surf_s = surf_def_v(l)%start_index(j,i)
2713	surf_e = surf_def_v(l)%end_index(j,i)
2714	DO m = surf_s, surf_e
2715	k = surf_def_v(l)%k(m)
2716	vsus = surf_def_v(l)%mom_flux_tke(0,m)
2717	wsus = surf_def_v(l)%mom_flux_tke(1,m)
2718
2719	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2720	* 0.5_wp * dx
2721	!
2722	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2723	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2724	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2725	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2726	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2727	ENDDO
2728	!
2729	!-- Natural surfaces
2730	surf_s = surf_lsm_v(l)%start_index(j,i)
2731	surf_e = surf_lsm_v(l)%end_index(j,i)
2732	DO m = surf_s, surf_e
2733	k = surf_lsm_v(l)%k(m)
2734	vsus = surf_lsm_v(l)%mom_flux_tke(0,m)
2735	wsus = surf_lsm_v(l)%mom_flux_tke(1,m)
2736
2737	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2738	* 0.5_wp * dx
2739	!
2740	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2741	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2742	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2743	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2744	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2745	ENDDO
2746	!
2747	!-- Urban surfaces
2748	surf_s = surf_usm_v(l)%start_index(j,i)
2749	surf_e = surf_usm_v(l)%end_index(j,i)
2750	DO m = surf_s, surf_e
2751	k = surf_usm_v(l)%k(m)
2752	vsus = surf_usm_v(l)%mom_flux_tke(0,m)
2753	wsus = surf_usm_v(l)%mom_flux_tke(1,m)
2754
2755	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2756	* 0.5_wp * dx
2757	!
2758	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2759	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2760	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2761	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2762	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2763	ENDDO
2764	ENDDO
2765	!
2766	!-- Compute gradients at upward-facing surfaces
2767	surf_s = surf_def_h(0)%start_index(j,i)
2768	surf_e = surf_def_h(0)%end_index(j,i)
2769	DO m = surf_s, surf_e
2770	k = surf_def_h(0)%k(m)
2771	!
2772	!-- Please note, actually, an interpolation of u_0 and v_0
2773	!-- onto the grid center would be required. However, this
2774	!-- would require several data transfers between 2D-grid and
2775	!-- wall type. The effect of this missing interpolation is
2776	!-- negligible. (See also production_e_init).
2777	dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
2778	dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
2779
2780	ENDDO
2781	!
2782	!-- Natural surfaces
2783	surf_s = surf_lsm_h%start_index(j,i)
2784	surf_e = surf_lsm_h%end_index(j,i)
2785	DO m = surf_s, surf_e
2786	k = surf_lsm_h%k(m)
2787
2788	dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
2789	dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
2790
2791	ENDDO
2792	!
2793	!-- Urban surfaces
2794	surf_s = surf_usm_h%start_index(j,i)
2795	surf_e = surf_usm_h%end_index(j,i)
2796	DO m = surf_s, surf_e
2797	k = surf_usm_h%k(m)
2798
2799	dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
2800	dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
2801
2802	ENDDO
2803	!
2804	!-- Compute gradients at downward-facing walls, only for
2805	!-- non-natural default surfaces
2806	surf_s = surf_def_h(1)%start_index(j,i)
2807	surf_e = surf_def_h(1)%end_index(j,i)
2808	DO m = surf_s, surf_e
2809	k = surf_def_h(1)%k(m)
2810
2811	dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
2812	dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
2813
2814	ENDDO
2815
2816
2817	ENDIF
2818
2819
2820	DO k = nzb+1, nzt
2821
2822	def = 2.0_wp * ( dudx(k)2 + dvdy(k)2 + dwdz(k)**2 ) + &
2823	dudy(k)2 + dvdx(k)2 + dwdx(k)**2 + &
2824	dwdy(k)2 + dudz(k)2 + dvdz(k)**2 + &
2825	2.0_wp * ( dvdx(k)dudy(k) + dwdx(k)dudz(k) + &
2826	dwdy(k)*dvdz(k) )
2827
2828	IF ( def < 0.0_wp ) def = 0.0_wp
2829
2830	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
2831
2832	IF ( .NOT. diss_production ) THEN
2833
2834	!-- Compute temdency for TKE-production from shear
2835	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
2836
2837	ELSE
2838
2839	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
2840	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * &
2841	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
2842
2843	ENDIF
2844
2845	ENDDO
2846
2847
2848	ENDDO
2849	ENDDO
2850
2851	!
2852	!-- If required, calculate TKE production by buoyancy
2853	IF ( .NOT. neutral ) THEN
2854
2855	IF ( .NOT. humidity ) THEN
2856
2857	IF ( ocean_mode ) THEN
2858	!
2859	!-- So far in the ocean no special treatment of density flux
2860	!-- in the bottom and top surface layer
2861	DO i = nxl, nxr
2862	DO j = nys, nyn
2863
2864	DO k = nzb+1, nzt
2865	tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
2866	ENDDO
2867	!
2868	!-- Treatment of near-surface grid points, at up- and down-
2869	!-- ward facing surfaces
2870	IF ( use_surface_fluxes ) THEN
2871	DO l = 0, 1
2872	surf_s = surf_def_h(l)%start_index(j,i)
2873	surf_e = surf_def_h(l)%end_index(j,i)
2874	DO m = surf_s, surf_e
2875	k = surf_def_h(l)%k(m)
2876	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
2877	ENDDO
2878	ENDDO
2879	ENDIF
2880
2881	IF ( use_top_fluxes ) THEN
2882	surf_s = surf_def_h(2)%start_index(j,i)
2883	surf_e = surf_def_h(2)%end_index(j,i)
2884	DO m = surf_s, surf_e
2885	k = surf_def_h(2)%k(m)
2886	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
2887	ENDDO
2888	ENDIF
2889
2890	IF ( .NOT. diss_production ) THEN
2891
2892	!-- Compute temdency for TKE-production from shear
2893	DO k = nzb+1, nzt
2894	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2895	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2896	MERGE( rho_reference, prho(k,j,i), &
2897	use_single_reference_value ) )
2898	ENDDO
2899
2900	ELSE
2901
2902	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
2903	DO k = nzb+1, nzt
2904	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2905	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2906	MERGE( rho_reference, prho(k,j,i), &
2907	use_single_reference_value ) ) * &
2908	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
2909	c_3
2910	ENDDO
2911
2912	ENDIF
2913
2914	ENDDO
2915	ENDDO
2916
2917	ELSE ! or IF ( .NOT. ocean_mode ) THEN
2918
2919	DO i = nxl, nxr
2920	DO j = nys, nyn
2921
2922	DO k = nzb+1, nzt
2923	tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
2924	ENDDO
2925
2926	IF ( use_surface_fluxes ) THEN
2927	!
2928	!-- Default surfaces, up- and downward-facing
2929	DO l = 0, 1
2930	surf_s = surf_def_h(l)%start_index(j,i)
2931	surf_e = surf_def_h(l)%end_index(j,i)
2932	DO m = surf_s, surf_e
2933	k = surf_def_h(l)%k(m)
2934	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
2935	ENDDO
2936	ENDDO
2937	!
2938	!-- Natural surfaces
2939	surf_s = surf_lsm_h%start_index(j,i)
2940	surf_e = surf_lsm_h%end_index(j,i)
2941	DO m = surf_s, surf_e
2942	k = surf_lsm_h%k(m)
2943	tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
2944	ENDDO
2945	!
2946	!-- Urban surfaces
2947	surf_s = surf_usm_h%start_index(j,i)
2948	surf_e = surf_usm_h%end_index(j,i)
2949	DO m = surf_s, surf_e
2950	k = surf_usm_h%k(m)
2951	tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
2952	ENDDO
2953	ENDIF
2954
2955	IF ( use_top_fluxes ) THEN
2956	surf_s = surf_def_h(2)%start_index(j,i)
2957	surf_e = surf_def_h(2)%end_index(j,i)
2958	DO m = surf_s, surf_e
2959	k = surf_def_h(2)%k(m)
2960	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
2961	ENDDO
2962	ENDIF
2963
2964	IF ( .NOT. diss_production ) THEN
2965
2966	!-- Compute temdency for TKE-production from shear
2967	DO k = nzb+1, nzt
2968	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2969	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2970	MERGE( pt_reference, pt(k,j,i), &
2971	use_single_reference_value ) )
2972	ENDDO
2973
2974	ELSE
2975
2976	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
2977	DO k = nzb+1, nzt
2978	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2979	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2980	MERGE( pt_reference, pt(k,j,i), &
2981	use_single_reference_value ) ) * &
2982	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
2983	c_3
2984	ENDDO
2985
2986	ENDIF
2987
2988	ENDDO
2989	ENDDO
2990
2991	ENDIF ! from IF ( .NOT. ocean_mode )
2992
2993	ELSE ! or IF ( humidity ) THEN
2994
2995	DO i = nxl, nxr
2996	DO j = nys, nyn
2997
2998	DO k = nzb+1, nzt
2999
3000	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3001	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3002	k2 = 0.61_wp * pt(k,j,i)
3003	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3004	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3005	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
3006	) * dd2zu(k)
3007	ELSE IF ( bulk_cloud_model ) THEN
3008	IF ( ql(k,j,i) == 0.0_wp ) THEN
3009	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3010	k2 = 0.61_wp * pt(k,j,i)
3011	ELSE
3012	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3013	temp = theta * exner(k)
3014	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3015	( q(k,j,i) - ql(k,j,i) ) * &
3016	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3017	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3018	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3019	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3020	ENDIF
3021	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3022	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3023	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
3024	) * dd2zu(k)
3025	ELSE IF ( cloud_droplets ) THEN
3026	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3027	k2 = 0.61_wp * pt(k,j,i)
3028	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3029	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3030	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
3031	pt(k,j,i) * ( ql(k+1,j,i) - &
3032	ql(k-1,j,i) ) ) * dd2zu(k)
3033	ENDIF
3034
3035	ENDDO
3036
3037	IF ( use_surface_fluxes ) THEN
3038	!
3039	!-- Treat horizontal default surfaces
3040	DO l = 0, 1
3041	surf_s = surf_def_h(l)%start_index(j,i)
3042	surf_e = surf_def_h(l)%end_index(j,i)
3043	DO m = surf_s, surf_e
3044	k = surf_def_h(l)%k(m)
3045
3046	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3047	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3048	k2 = 0.61_wp * pt(k,j,i)
3049	ELSE IF ( bulk_cloud_model ) THEN
3050	IF ( ql(k,j,i) == 0.0_wp ) THEN
3051	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3052	k2 = 0.61_wp * pt(k,j,i)
3053	ELSE
3054	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3055	temp = theta * exner(k)
3056	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3057	( q(k,j,i) - ql(k,j,i) ) * &
3058	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3059	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3060	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3061	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3062	ENDIF
3063	ELSE IF ( cloud_droplets ) THEN
3064	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3065	k2 = 0.61_wp * pt(k,j,i)
3066	ENDIF
3067
3068	tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + &
3069	k2 * surf_def_h(l)%qsws(m) &
3070	) * drho_air_zw(k-1)
3071	ENDDO
3072	ENDDO
3073	!
3074	!-- Treat horizontal natural surfaces
3075	surf_s = surf_lsm_h%start_index(j,i)
3076	surf_e = surf_lsm_h%end_index(j,i)
3077	DO m = surf_s, surf_e
3078	k = surf_lsm_h%k(m)
3079
3080	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3081	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3082	k2 = 0.61_wp * pt(k,j,i)
3083	ELSE IF ( bulk_cloud_model ) THEN
3084	IF ( ql(k,j,i) == 0.0_wp ) THEN
3085	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3086	k2 = 0.61_wp * pt(k,j,i)
3087	ELSE
3088	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3089	temp = theta * exner(k)
3090	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3091	( q(k,j,i) - ql(k,j,i) ) * &
3092	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3093	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3094	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3095	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3096	ENDIF
3097	ELSE IF ( cloud_droplets ) THEN
3098	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3099	k2 = 0.61_wp * pt(k,j,i)
3100	ENDIF
3101
3102	tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + &
3103	k2 * surf_lsm_h%qsws(m) &
3104	) * drho_air_zw(k-1)
3105	ENDDO
3106	!
3107	!-- Treat horizontal urban surfaces
3108	surf_s = surf_usm_h%start_index(j,i)
3109	surf_e = surf_usm_h%end_index(j,i)
3110	DO m = surf_s, surf_e
3111	k = surf_usm_h%k(m)
3112
3113	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3114	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3115	k2 = 0.61_wp * pt(k,j,i)
3116	ELSE IF ( bulk_cloud_model ) THEN
3117	IF ( ql(k,j,i) == 0.0_wp ) THEN
3118	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3119	k2 = 0.61_wp * pt(k,j,i)
3120	ELSE
3121	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3122	temp = theta * exner(k)
3123	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3124	( q(k,j,i) - ql(k,j,i) ) * &
3125	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3126	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3127	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3128	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3129	ENDIF
3130	ELSE IF ( cloud_droplets ) THEN
3131	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3132	k2 = 0.61_wp * pt(k,j,i)
3133	ENDIF
3134
3135	tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + &
3136	k2 * surf_usm_h%qsws(m) &
3137	) * drho_air_zw(k-1)
3138	ENDDO
3139
3140	ENDIF ! from IF ( use_surface_fluxes ) THEN
3141
3142	IF ( use_top_fluxes ) THEN
3143
3144	surf_s = surf_def_h(2)%start_index(j,i)
3145	surf_e = surf_def_h(2)%end_index(j,i)
3146	DO m = surf_s, surf_e
3147	k = surf_def_h(2)%k(m)
3148
3149	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3150	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3151	k2 = 0.61_wp * pt(k,j,i)
3152	ELSE IF ( bulk_cloud_model ) THEN
3153	IF ( ql(k,j,i) == 0.0_wp ) THEN
3154	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3155	k2 = 0.61_wp * pt(k,j,i)
3156	ELSE
3157	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3158	temp = theta * exner(k)
3159	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3160	( q(k,j,i) - ql(k,j,i) ) * &
3161	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3162	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3163	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3164	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3165	ENDIF
3166	ELSE IF ( cloud_droplets ) THEN
3167	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3168	k2 = 0.61_wp * pt(k,j,i)
3169	ENDIF
3170
3171	tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + &
3172	k2 * surf_def_h(2)%qsws(m) &
3173	) * drho_air_zw(k)
3174
3175	ENDDO
3176
3177	ENDIF ! from IF ( use_top_fluxes ) THEN
3178
3179	IF ( .NOT. diss_production ) THEN
3180
3181	!-- Compute temdency for TKE-production from shear
3182	DO k = nzb+1, nzt
3183	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3184	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3185	MERGE( vpt_reference, vpt(k,j,i), &
3186	use_single_reference_value ) )
3187	ENDDO
3188
3189	ELSE
3190
3191	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
3192	DO k = nzb+1, nzt
3193	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3194	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3195	MERGE( vpt_reference, vpt(k,j,i), &
3196	use_single_reference_value ) ) * &
3197	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3198	c_3
3199	ENDDO
3200
3201	ENDIF
3202
3203	ENDDO
3204	ENDDO
3205
3206	ENDIF
3207
3208	ENDIF
3209
3210	END SUBROUTINE production_e
3211
3212
3213	!------------------------------------------------------------------------------!
3214	! Description:
3215	! ------------
3216	!> Production terms (shear + buoyancy) of the TKE.
3217	!> Cache-optimized version
3218	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
3219	!> not considered well!
3220	!> @todo non-neutral case is not yet considered for RANS mode
3221	!------------------------------------------------------------------------------!
3222	SUBROUTINE production_e_ij( i, j, diss_production )
3223
3224	USE arrays_3d, &
3225	ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
3226
3227	USE control_parameters, &
3228	ONLY: cloud_droplets, constant_flux_layer, neutral, &
3229	rho_reference, use_single_reference_value, use_surface_fluxes, &
3230	use_top_fluxes
3231
3232	USE grid_variables, &
3233	ONLY: ddx, dx, ddy, dy
3234
3235	USE bulk_cloud_model_mod, &
3236	ONLY: bulk_cloud_model
3237
3238	USE surface_mod, &
3239	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
3240	surf_usm_v
3241
3242	IMPLICIT NONE
3243
3244	LOGICAL :: diss_production
3245
3246	INTEGER(iwp) :: i !< running index x-direction
3247	INTEGER(iwp) :: j !< running index y-direction
3248	INTEGER(iwp) :: k !< running index z-direction
3249	INTEGER(iwp) :: l !< running index for different surface type orientation
3250	INTEGER(iwp) :: m !< running index surface elements
3251	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
3252	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
3253	INTEGER(iwp) :: flag_nr !< number of masking flag
3254
3255	REAL(wp) :: def !<
3256	REAL(wp) :: flag !< flag to mask topography
3257	REAL(wp) :: k1 !<
3258	REAL(wp) :: k2 !<
3259	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
3260	REAL(wp) :: theta !<
3261	REAL(wp) :: temp !<
3262	REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation
3263	REAL(wp) :: usvs !< momentum flux u"v"
3264	REAL(wp) :: vsus !< momentum flux v"u"
3265	REAL(wp) :: wsus !< momentum flux w"u"
3266	REAL(wp) :: wsvs !< momentum flux w"v"
3267
3268	REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction
3269	REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction
3270	REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction
3271	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction
3272	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction
3273	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction
3274	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction
3275	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction
3276	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction
3277	REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary flux-array in z-direction
3278
3279
3280
3281	!
3282	!-- Calculate TKE production by shear. Calculate gradients at all grid
3283	!-- points first, gradients at surface-bounded grid points will be
3284	!-- overwritten further below.
3285	DO k = nzb+1, nzt
3286
3287	dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
3288	dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
3289	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
3290	dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
3291	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
3292
3293	dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
3294	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
3295	dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
3296	dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
3297	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
3298
3299	dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
3300	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
3301	dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
3302	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
3303	dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
3304
3305	ENDDO
3306
3307	flag_nr = 29
3308
3309	IF ( constant_flux_layer ) THEN
3310
3311	flag_nr = 0
3312
3313	!-- Position beneath wall
3314	!-- (2) - Will allways be executed.
3315	!-- 'bottom and wall: use u_0,v_0 and wall functions'
3316	!
3317	!-- Compute gradients at north- and south-facing surfaces.
3318	!-- First, for default surfaces, then for urban surfaces.
3319	!-- Note, so far no natural vertical surfaces implemented
3320	DO l = 0, 1
3321	surf_s = surf_def_v(l)%start_index(j,i)
3322	surf_e = surf_def_v(l)%end_index(j,i)
3323	DO m = surf_s, surf_e
3324	k = surf_def_v(l)%k(m)
3325	usvs = surf_def_v(l)%mom_flux_tke(0,m)
3326	wsvs = surf_def_v(l)%mom_flux_tke(1,m)
3327
3328	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3329	* 0.5_wp * dy
3330	!
3331	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3332	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3333	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3334	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3335	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3336	ENDDO
3337	!
3338	!-- Natural surfaces
3339	surf_s = surf_lsm_v(l)%start_index(j,i)
3340	surf_e = surf_lsm_v(l)%end_index(j,i)
3341	DO m = surf_s, surf_e
3342	k = surf_lsm_v(l)%k(m)
3343	usvs = surf_lsm_v(l)%mom_flux_tke(0,m)
3344	wsvs = surf_lsm_v(l)%mom_flux_tke(1,m)
3345
3346	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3347	* 0.5_wp * dy
3348	!
3349	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3350	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3351	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3352	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3353	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3354	ENDDO
3355	!
3356	!-- Urban surfaces
3357	surf_s = surf_usm_v(l)%start_index(j,i)
3358	surf_e = surf_usm_v(l)%end_index(j,i)
3359	DO m = surf_s, surf_e
3360	k = surf_usm_v(l)%k(m)
3361	usvs = surf_usm_v(l)%mom_flux_tke(0,m)
3362	wsvs = surf_usm_v(l)%mom_flux_tke(1,m)
3363
3364	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3365	* 0.5_wp * dy
3366	!
3367	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3368	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3369	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3370	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3371	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3372	ENDDO
3373	ENDDO
3374	!
3375	!-- Compute gradients at east- and west-facing walls
3376	DO l = 2, 3
3377	surf_s = surf_def_v(l)%start_index(j,i)
3378	surf_e = surf_def_v(l)%end_index(j,i)
3379	DO m = surf_s, surf_e
3380	k = surf_def_v(l)%k(m)
3381	vsus = surf_def_v(l)%mom_flux_tke(0,m)
3382	wsus = surf_def_v(l)%mom_flux_tke(1,m)
3383
3384	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3385	* 0.5_wp * dx
3386	!
3387	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3388	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3389	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3390	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3391	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3392	ENDDO
3393	!
3394	!-- Natural surfaces
3395	surf_s = surf_lsm_v(l)%start_index(j,i)
3396	surf_e = surf_lsm_v(l)%end_index(j,i)
3397	DO m = surf_s, surf_e
3398	k = surf_lsm_v(l)%k(m)
3399	vsus = surf_lsm_v(l)%mom_flux_tke(0,m)
3400	wsus = surf_lsm_v(l)%mom_flux_tke(1,m)
3401
3402	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3403	* 0.5_wp * dx
3404	!
3405	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3406	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3407	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3408	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3409	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3410	ENDDO
3411	!
3412	!-- Urban surfaces
3413	surf_s = surf_usm_v(l)%start_index(j,i)
3414	surf_e = surf_usm_v(l)%end_index(j,i)
3415	DO m = surf_s, surf_e
3416	k = surf_usm_v(l)%k(m)
3417	vsus = surf_usm_v(l)%mom_flux_tke(0,m)
3418	wsus = surf_usm_v(l)%mom_flux_tke(1,m)
3419
3420	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3421	* 0.5_wp * dx
3422	!
3423	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3424	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3425	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3426	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3427	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3428	ENDDO
3429	ENDDO
3430	!
3431	!-- Compute gradients at upward-facing surfaces
3432	surf_s = surf_def_h(0)%start_index(j,i)
3433	surf_e = surf_def_h(0)%end_index(j,i)
3434	DO m = surf_s, surf_e
3435	k = surf_def_h(0)%k(m)
3436	!
3437	!-- Please note, actually, an interpolation of u_0 and v_0
3438	!-- onto the grid center would be required. However, this
3439	!-- would require several data transfers between 2D-grid and
3440	!-- wall type. The effect of this missing interpolation is
3441	!-- negligible. (See also production_e_init).
3442	dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
3443	dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
3444
3445	ENDDO
3446	!
3447	!-- Natural surfaces
3448	surf_s = surf_lsm_h%start_index(j,i)
3449	surf_e = surf_lsm_h%end_index(j,i)
3450	DO m = surf_s, surf_e
3451	k = surf_lsm_h%k(m)
3452
3453	dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
3454	dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
3455
3456	ENDDO
3457	!
3458	!-- Urban surfaces
3459	surf_s = surf_usm_h%start_index(j,i)
3460	surf_e = surf_usm_h%end_index(j,i)
3461	DO m = surf_s, surf_e
3462	k = surf_usm_h%k(m)
3463
3464	dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
3465	dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
3466
3467	ENDDO
3468	!
3469	!-- Compute gradients at downward-facing walls, only for
3470	!-- non-natural default surfaces
3471	surf_s = surf_def_h(1)%start_index(j,i)
3472	surf_e = surf_def_h(1)%end_index(j,i)
3473	DO m = surf_s, surf_e
3474	k = surf_def_h(1)%k(m)
3475
3476	dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
3477	dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
3478
3479	ENDDO
3480
3481	ENDIF
3482
3483	DO k = nzb+1, nzt
3484
3485	def = 2.0_wp * ( dudx(k)2 + dvdy(k)2 + dwdz(k)**2 ) + &
3486	dudy(k)2 + dvdx(k)2 + dwdx(k)**2 + &
3487	dwdy(k)2 + dudz(k)2 + dvdz(k)**2 + &
3488	2.0_wp * ( dvdx(k)dudy(k) + dwdx(k)dudz(k) + &
3489	dwdy(k)*dvdz(k) )
3490
3491	IF ( def < 0.0_wp ) def = 0.0_wp
3492
3493	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
3494
3495	IF ( .NOT. diss_production ) THEN
3496
3497	!-- Compute temdency for TKE-production from shear
3498	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
3499
3500	ELSE
3501
3502	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
3503	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * &
3504	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
3505
3506	ENDIF
3507
3508	ENDDO
3509
3510
3511	!
3512	!-- If required, calculate TKE production by buoyancy
3513	IF ( .NOT. neutral ) THEN
3514
3515	IF ( .NOT. humidity ) THEN
3516
3517	IF ( ocean_mode ) THEN
3518	!
3519	!-- So far in the ocean no special treatment of density flux
3520	!-- in the bottom and top surface layer
3521
3522	DO k = nzb+1, nzt
3523	tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
3524	ENDDO
3525	!
3526	!-- Treatment of near-surface grid points, at up- and down-
3527	!-- ward facing surfaces
3528	IF ( use_surface_fluxes ) THEN
3529	DO l = 0, 1
3530	surf_s = surf_def_h(l)%start_index(j,i)
3531	surf_e = surf_def_h(l)%end_index(j,i)
3532	DO m = surf_s, surf_e
3533	k = surf_def_h(l)%k(m)
3534	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3535	ENDDO
3536	ENDDO
3537	ENDIF
3538
3539	IF ( use_top_fluxes ) THEN
3540	surf_s = surf_def_h(2)%start_index(j,i)
3541	surf_e = surf_def_h(2)%end_index(j,i)
3542	DO m = surf_s, surf_e
3543	k = surf_def_h(2)%k(m)
3544	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3545	ENDDO
3546	ENDIF
3547
3548	IF ( .NOT. diss_production ) THEN
3549
3550	!-- Compute temdency for TKE-production from shear
3551	DO k = nzb+1, nzt
3552	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3553	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3554	MERGE( rho_reference, prho(k,j,i), &
3555	use_single_reference_value ) )
3556	ENDDO
3557
3558	ELSE
3559
3560	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
3561	DO k = nzb+1, nzt
3562	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3563	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3564	MERGE( rho_reference, prho(k,j,i), &
3565	use_single_reference_value ) ) * &
3566	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3567	c_3
3568	ENDDO
3569
3570	ENDIF
3571
3572
3573	ELSE ! or IF ( .NOT. ocean_mode ) THEN
3574
3575	DO k = nzb+1, nzt
3576	tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
3577	ENDDO
3578
3579	IF ( use_surface_fluxes ) THEN
3580	!
3581	!-- Default surfaces, up- and downward-facing
3582	DO l = 0, 1
3583	surf_s = surf_def_h(l)%start_index(j,i)
3584	surf_e = surf_def_h(l)%end_index(j,i)
3585	DO m = surf_s, surf_e
3586	k = surf_def_h(l)%k(m)
3587	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3588	ENDDO
3589	ENDDO
3590	!
3591	!-- Natural surfaces
3592	surf_s = surf_lsm_h%start_index(j,i)
3593	surf_e = surf_lsm_h%end_index(j,i)
3594	DO m = surf_s, surf_e
3595	k = surf_lsm_h%k(m)
3596	tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
3597	ENDDO
3598	!
3599	!-- Urban surfaces
3600	surf_s = surf_usm_h%start_index(j,i)
3601	surf_e = surf_usm_h%end_index(j,i)
3602	DO m = surf_s, surf_e
3603	k = surf_usm_h%k(m)
3604	tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
3605	ENDDO
3606	ENDIF
3607
3608	IF ( use_top_fluxes ) THEN
3609	surf_s = surf_def_h(2)%start_index(j,i)
3610	surf_e = surf_def_h(2)%end_index(j,i)
3611	DO m = surf_s, surf_e
3612	k = surf_def_h(2)%k(m)
3613	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3614	ENDDO
3615	ENDIF
3616
3617	IF ( .NOT. diss_production ) THEN
3618
3619	!-- Compute temdency for TKE-production from shear
3620	DO k = nzb+1, nzt
3621	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3622	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3623	MERGE( pt_reference, pt(k,j,i), &
3624	use_single_reference_value ) )
3625	ENDDO
3626
3627	ELSE
3628
3629	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
3630	DO k = nzb+1, nzt
3631	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3632	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3633	MERGE( pt_reference, pt(k,j,i), &
3634	use_single_reference_value ) ) * &
3635	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3636	c_3
3637	ENDDO
3638
3639	ENDIF
3640
3641	ENDIF ! from IF ( .NOT. ocean_mode )
3642
3643	ELSE ! or IF ( humidity ) THEN
3644
3645	DO k = nzb+1, nzt
3646
3647	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3648	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3649	k2 = 0.61_wp * pt(k,j,i)
3650	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3651	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3652	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
3653	) * dd2zu(k)
3654	ELSE IF ( bulk_cloud_model ) THEN
3655	IF ( ql(k,j,i) == 0.0_wp ) THEN
3656	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3657	k2 = 0.61_wp * pt(k,j,i)
3658	ELSE
3659	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3660	temp = theta * exner(k)
3661	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3662	( q(k,j,i) - ql(k,j,i) ) * &
3663	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3664	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3665	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3666	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3667	ENDIF
3668	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3669	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3670	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
3671	) * dd2zu(k)
3672	ELSE IF ( cloud_droplets ) THEN
3673	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3674	k2 = 0.61_wp * pt(k,j,i)
3675	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3676	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3677	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
3678	pt(k,j,i) * ( ql(k+1,j,i) - &
3679	ql(k-1,j,i) ) ) * dd2zu(k)
3680	ENDIF
3681
3682	ENDDO
3683
3684	IF ( use_surface_fluxes ) THEN
3685	!
3686	!-- Treat horizontal default surfaces
3687	DO l = 0, 1
3688	surf_s = surf_def_h(l)%start_index(j,i)
3689	surf_e = surf_def_h(l)%end_index(j,i)
3690	DO m = surf_s, surf_e
3691	k = surf_def_h(l)%k(m)
3692
3693	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3694	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3695	k2 = 0.61_wp * pt(k,j,i)
3696	ELSE IF ( bulk_cloud_model ) THEN
3697	IF ( ql(k,j,i) == 0.0_wp ) THEN
3698	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3699	k2 = 0.61_wp * pt(k,j,i)
3700	ELSE
3701	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3702	temp = theta * exner(k)
3703	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3704	( q(k,j,i) - ql(k,j,i) ) * &
3705	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3706	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3707	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3708	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3709	ENDIF
3710	ELSE IF ( cloud_droplets ) THEN
3711	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3712	k2 = 0.61_wp * pt(k,j,i)
3713	ENDIF
3714
3715	tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + &
3716	k2 * surf_def_h(l)%qsws(m) &
3717	) * drho_air_zw(k-1)
3718	ENDDO
3719	ENDDO
3720	!
3721	!-- Treat horizontal natural surfaces
3722	surf_s = surf_lsm_h%start_index(j,i)
3723	surf_e = surf_lsm_h%end_index(j,i)
3724	DO m = surf_s, surf_e
3725	k = surf_lsm_h%k(m)
3726
3727	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3728	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3729	k2 = 0.61_wp * pt(k,j,i)
3730	ELSE IF ( bulk_cloud_model ) THEN
3731	IF ( ql(k,j,i) == 0.0_wp ) THEN
3732	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3733	k2 = 0.61_wp * pt(k,j,i)
3734	ELSE
3735	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3736	temp = theta * exner(k)
3737	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3738	( q(k,j,i) - ql(k,j,i) ) * &
3739	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3740	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3741	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3742	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3743	ENDIF
3744	ELSE IF ( cloud_droplets ) THEN
3745	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3746	k2 = 0.61_wp * pt(k,j,i)
3747	ENDIF
3748
3749	tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + &
3750	k2 * surf_lsm_h%qsws(m) &
3751	) * drho_air_zw(k-1)
3752	ENDDO
3753	!
3754	!-- Treat horizontal urban surfaces
3755	surf_s = surf_usm_h%start_index(j,i)
3756	surf_e = surf_usm_h%end_index(j,i)
3757	DO m = surf_s, surf_e
3758	k = surf_usm_h%k(m)
3759
3760	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3761	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3762	k2 = 0.61_wp * pt(k,j,i)
3763	ELSE IF ( bulk_cloud_model ) THEN
3764	IF ( ql(k,j,i) == 0.0_wp ) THEN
3765	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3766	k2 = 0.61_wp * pt(k,j,i)
3767	ELSE
3768	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3769	temp = theta * exner(k)
3770	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3771	( q(k,j,i) - ql(k,j,i) ) * &
3772	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3773	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3774	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3775	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3776	ENDIF
3777	ELSE IF ( cloud_droplets ) THEN
3778	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3779	k2 = 0.61_wp * pt(k,j,i)
3780	ENDIF
3781
3782	tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + &
3783	k2 * surf_usm_h%qsws(m) &
3784	) * drho_air_zw(k-1)
3785	ENDDO
3786
3787	ENDIF ! from IF ( use_surface_fluxes ) THEN
3788
3789	IF ( use_top_fluxes ) THEN
3790
3791	surf_s = surf_def_h(2)%start_index(j,i)
3792	surf_e = surf_def_h(2)%end_index(j,i)
3793	DO m = surf_s, surf_e
3794	k = surf_def_h(2)%k(m)
3795
3796	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3797	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3798	k2 = 0.61_wp * pt(k,j,i)
3799	ELSE IF ( bulk_cloud_model ) THEN
3800	IF ( ql(k,j,i) == 0.0_wp ) THEN
3801	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3802	k2 = 0.61_wp * pt(k,j,i)
3803	ELSE
3804	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3805	temp = theta * exner(k)
3806	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3807	( q(k,j,i) - ql(k,j,i) ) * &
3808	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3809	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3810	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3811	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3812	ENDIF
3813	ELSE IF ( cloud_droplets ) THEN
3814	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3815	k2 = 0.61_wp * pt(k,j,i)
3816	ENDIF
3817
3818	tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + &
3819	k2 * surf_def_h(2)%qsws(m) &
3820	) * drho_air_zw(k)
3821
3822	ENDDO
3823
3824	ENDIF ! from IF ( use_top_fluxes ) THEN
3825
3826	IF ( .NOT. diss_production ) THEN
3827
3828	!-- Compute temdency for TKE-production from shear
3829	DO k = nzb+1, nzt
3830	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3831	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3832	MERGE( vpt_reference, vpt(k,j,i), &
3833	use_single_reference_value ) )
3834	ENDDO
3835
3836	ELSE
3837
3838	!-- RANS mode: Compute temdency for dissipation-rate-production from shear
3839	DO k = nzb+1, nzt
3840	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3841	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3842	MERGE( vpt_reference, vpt(k,j,i), &
3843	use_single_reference_value ) ) * &
3844	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3845	c_3
3846	ENDDO
3847
3848	ENDIF
3849
3850	ENDIF
3851
3852	ENDIF
3853
3854	END SUBROUTINE production_e_ij
3855
3856
3857	!------------------------------------------------------------------------------!
3858	! Description:
3859	! ------------
3860	!> Diffusion and dissipation terms for the TKE.
3861	!> Vector-optimized version
3862	!------------------------------------------------------------------------------!
3863	SUBROUTINE diffusion_e( var, var_reference )
3864
3865	USE arrays_3d, &
3866	ONLY: ddzu, ddzw, drho_air, rho_air_zw
3867
3868	USE grid_variables, &
3869	ONLY: ddx2, ddy2
3870
3871	USE bulk_cloud_model_mod, &
3872	ONLY: collision_turbulence
3873
3874	USE particle_attributes, &
3875	ONLY: use_sgs_for_particles, wang_kernel
3876
3877	USE surface_mod, &
3878	ONLY : bc_h
3879
3880	IMPLICIT NONE
3881
3882	INTEGER(iwp) :: i !< running index x direction
3883	INTEGER(iwp) :: j !< running index y direction
3884	INTEGER(iwp) :: k !< running index z direction
3885	INTEGER(iwp) :: m !< running index surface elements
3886
3887	REAL(wp) :: flag !< flag to mask topography
3888	REAL(wp) :: l !< mixing length
3889	REAL(wp) :: ll !< adjusted l
3890	REAL(wp) :: var_reference !< reference temperature
3891
3892	#if defined( __nopointer )
3893	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature
3894	#else
3895	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
3896	#endif
3897	REAL(wp), DIMENSION(nzb+1:nzt,nys:nyn) :: dissipation !< TKE dissipation
3898
3899
3900	!
3901	!-- Calculate the tendency terms
3902	DO i = nxl, nxr
3903	DO j = nys, nyn
3904	DO k = nzb+1, nzt
3905	!
3906	!-- Predetermine flag to mask topography
3907	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3908
3909	!
3910	!-- Calculate dissipation
3911	IF ( les_dynamic .OR. les_mw ) THEN
3912
3913	CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
3914
3915	dissipation(k,j) = ( 0.19_wp + 0.74_wp * l / ll ) &
3916	* e(k,j,i) * SQRT( e(k,j,i) ) / l
3917
3918	ELSEIF ( rans_tke_l ) THEN
3919
3920	CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
3921
3922	dissipation(k,j) = c_0*3 e(k,j,i) * SQRT( e(k,j,i) ) / ll
3923
3924	diss(k,j,i) = dissipation(k,j) * flag
3925
3926	ELSEIF ( rans_tke_e ) THEN
3927
3928	dissipation(k,j) = diss(k,j,i)
3929
3930	ENDIF
3931
3932	tend(k,j,i) = tend(k,j,i) + ( &
3933	( &
3934	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
3935	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
3936	) * ddx2 * flag &
3937	+ ( &
3938	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
3939	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
3940	) * ddy2 * flag &
3941	+ ( &
3942	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
3943	* rho_air_zw(k) &
3944	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
3945	* rho_air_zw(k-1) &
3946	) * ddzw(k) * drho_air(k) &
3947	) * flag * dsig_e &
3948	- dissipation(k,j) * flag
3949
3950	ENDDO
3951	ENDDO
3952
3953	!
3954	!-- Store dissipation if needed for calculating the sgs particle
3955	!-- velocities
3956	IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. &
3957	wang_kernel .OR. collision_turbulence ) ) THEN
3958	DO j = nys, nyn
3959	DO k = nzb+1, nzt
3960	diss(k,j,i) = dissipation(k,j) * MERGE( 1.0_wp, 0.0_wp, &
3961	BTEST( wall_flags_0(k,j,i), 0 ) )
3962	ENDDO
3963	ENDDO
3964	ENDIF
3965
3966	ENDDO
3967
3968	!
3969	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
3970	IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. &
3971	wang_kernel .OR. collision_turbulence ) ) THEN
3972	!
3973	!-- Upward facing surfaces
3974	DO m = 1, bc_h(0)%ns
3975	i = bc_h(0)%i(m)
3976	j = bc_h(0)%j(m)
3977	k = bc_h(0)%k(m)
3978	diss(k-1,j,i) = diss(k,j,i)
3979	ENDDO
3980	!
3981	!-- Downward facing surfaces
3982	DO m = 1, bc_h(1)%ns
3983	i = bc_h(1)%i(m)
3984	j = bc_h(1)%j(m)
3985	k = bc_h(1)%k(m)
3986	diss(k+1,j,i) = diss(k,j,i)
3987	ENDDO
3988
3989	ENDIF
3990
3991	END SUBROUTINE diffusion_e
3992
3993
3994	!------------------------------------------------------------------------------!
3995	! Description:
3996	! ------------
3997	!> Diffusion and dissipation terms for the TKE.
3998	!> Cache-optimized version
3999	!------------------------------------------------------------------------------!
4000	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
4001
4002	USE arrays_3d, &
4003	ONLY: ddzu, ddzw, drho_air, rho_air_zw
4004
4005	USE grid_variables, &
4006	ONLY: ddx2, ddy2
4007
4008	USE bulk_cloud_model_mod, &
4009	ONLY: collision_turbulence
4010
4011	USE particle_attributes, &
4012	ONLY: use_sgs_for_particles, wang_kernel
4013
4014	USE surface_mod, &
4015	ONLY : bc_h
4016
4017	IMPLICIT NONE
4018
4019	INTEGER(iwp) :: i !< running index x direction
4020	INTEGER(iwp) :: j !< running index y direction
4021	INTEGER(iwp) :: k !< running index z direction
4022	INTEGER(iwp) :: m !< running index surface elements
4023	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
4024	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
4025
4026	REAL(wp) :: flag !< flag to mask topography
4027	REAL(wp) :: l !< mixing length
4028	REAL(wp) :: ll !< adjusted l
4029	REAL(wp) :: var_reference !< reference temperature
4030
4031	#if defined( __nopointer )
4032	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature
4033	#else
4034	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4035	#endif
4036	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !< dissipation of TKE
4037
4038	!
4039	!-- Calculate the mixing length (for dissipation)
4040	DO k = nzb+1, nzt
4041	!
4042	!-- Predetermine flag to mask topography
4043	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4044
4045	!
4046	!-- Calculate dissipation...
4047	!-- ...in case of LES
4048	IF ( les_dynamic .OR. les_mw ) THEN
4049
4050	CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
4051
4052	dissipation(k) = ( 0.19_wp + 0.74_wp * l / ll ) &
4053	* e(k,j,i) * SQRT( e(k,j,i) ) / l
4054
4055	!
4056	!-- ...in case of RANS
4057	ELSEIF ( rans_tke_l ) THEN
4058
4059	CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
4060
4061	dissipation(k) = c_0*3 e(k,j,i) * SQRT( e(k,j,i) ) / ll
4062
4063	diss(k,j,i) = dissipation(k) * flag
4064
4065	ELSEIF ( rans_tke_e ) THEN
4066
4067	dissipation(k) = diss(k,j,i)
4068
4069	ENDIF
4070
4071	!
4072	!-- Calculate the tendency term
4073	tend(k,j,i) = tend(k,j,i) + ( &
4074	( &
4075	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
4076	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
4077	) * ddx2 &
4078	+ ( &
4079	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
4080	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
4081	) * ddy2 &
4082	+ ( &
4083	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
4084	* rho_air_zw(k) &
4085	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
4086	* rho_air_zw(k-1) &
4087	) * ddzw(k) * drho_air(k) &
4088	) * flag * dsig_e &
4089	- dissipation(k) * flag
4090
4091	ENDDO
4092
4093	!
4094	!-- Store dissipation if needed for calculating the sgs particle velocities
4095	IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. wang_kernel &
4096	.OR. collision_turbulence ) ) THEN
4097	DO k = nzb+1, nzt
4098	diss(k,j,i) = dissipation(k) * MERGE( 1.0_wp, 0.0_wp, &
4099	BTEST( wall_flags_0(k,j,i), 0 ) )
4100	ENDDO
4101	!
4102	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
4103	!-- For each surface type determine start and end index (in case of elevated
4104	!-- topography several up/downward facing surfaces may exist.
4105	surf_s = bc_h(0)%start_index(j,i)
4106	surf_e = bc_h(0)%end_index(j,i)
4107	DO m = surf_s, surf_e
4108	k = bc_h(0)%k(m)
4109	diss(k-1,j,i) = diss(k,j,i)
4110	ENDDO
4111	!
4112	!-- Downward facing surfaces
4113	surf_s = bc_h(1)%start_index(j,i)
4114	surf_e = bc_h(1)%end_index(j,i)
4115	DO m = surf_s, surf_e
4116	k = bc_h(1)%k(m)
4117	diss(k+1,j,i) = diss(k,j,i)
4118	ENDDO
4119	ENDIF
4120
4121	END SUBROUTINE diffusion_e_ij
4122
4123
4124	!------------------------------------------------------------------------------!
4125	! Description:
4126	! ------------
4127	!> Diffusion term for the TKE dissipation rate
4128	!> Vector-optimized version
4129	!------------------------------------------------------------------------------!
4130	SUBROUTINE diffusion_diss()
4131	USE arrays_3d, &
4132	ONLY: ddzu, ddzw, drho_air, rho_air_zw
4133
4134	USE grid_variables, &
4135	ONLY: ddx2, ddy2
4136
4137	IMPLICIT NONE
4138
4139	INTEGER(iwp) :: i !< running index x direction
4140	INTEGER(iwp) :: j !< running index y direction
4141	INTEGER(iwp) :: k !< running index z direction
4142
4143	REAL(wp) :: flag !< flag to mask topography
4144
4145	!
4146	!-- Calculate the tendency terms
4147	DO i = nxl, nxr
4148	DO j = nys, nyn
4149	DO k = nzb+1, nzt
4150
4151	!
4152	!-- Predetermine flag to mask topography
4153	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4154
4155	tend(k,j,i) = tend(k,j,i) + &
4156	( ( &
4157	( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) &
4158	- ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) &
4159	) * ddx2 &
4160	+ ( &
4161	( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) &
4162	- ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) &
4163	) * ddy2 &
4164	+ ( &
4165	( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) &
4166	* rho_air_zw(k) &
4167	- ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) &
4168	* rho_air_zw(k-1) &
4169	) * ddzw(k) * drho_air(k) &
4170	) * flag * dsig_diss &
4171	- c_2 * diss(k,j,i)**2 &
4172	/ ( e(k,j,i) + 1.0E-20_wp ) * flag
4173
4174	ENDDO
4175	ENDDO
4176	ENDDO
4177
4178	END SUBROUTINE diffusion_diss
4179
4180
4181	!------------------------------------------------------------------------------!
4182	! Description:
4183	! ------------
4184	!> Diffusion term for the TKE dissipation rate
4185	!> Cache-optimized version
4186	!------------------------------------------------------------------------------!
4187	SUBROUTINE diffusion_diss_ij( i, j )
4188
4189	USE arrays_3d, &
4190	ONLY: ddzu, ddzw, drho_air, rho_air_zw
4191
4192	USE grid_variables, &
4193	ONLY: ddx2, ddy2
4194
4195	IMPLICIT NONE
4196
4197	INTEGER(iwp) :: i !< running index x direction
4198	INTEGER(iwp) :: j !< running index y direction
4199	INTEGER(iwp) :: k !< running index z direction
4200
4201	REAL(wp) :: flag !< flag to mask topography
4202
4203	!
4204	!-- Calculate the mixing length (for dissipation)
4205	DO k = nzb+1, nzt
4206
4207	!
4208	!-- Predetermine flag to mask topography
4209	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4210
4211	!
4212	!-- Calculate the tendency term
4213	tend(k,j,i) = tend(k,j,i) + &
4214	( ( &
4215	( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) &
4216	- ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) &
4217	) * ddx2 &
4218	+ ( &
4219	( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) &
4220	- ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) &
4221	) * ddy2 &
4222	+ ( &
4223	( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) &
4224	* rho_air_zw(k) &
4225	- ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) &
4226	* rho_air_zw(k-1) &
4227	) * ddzw(k) * drho_air(k) &
4228	) * flag * dsig_diss &
4229	- c_2 * diss(k,j,i)*2 / ( e(k,j,i) + 1.0E-20_wp ) flag
4230
4231	ENDDO
4232
4233	END SUBROUTINE diffusion_diss_ij
4234
4235
4236	!------------------------------------------------------------------------------!
4237	! Description:
4238	! ------------
4239	!> Calculate mixing length for LES mode.
4240	!------------------------------------------------------------------------------!
4241	SUBROUTINE mixing_length_les( i, j, k, l, ll, var, var_reference )
4242
4243	USE arrays_3d, &
4244	ONLY: dd2zu
4245
4246	USE control_parameters, &
4247	ONLY: atmos_ocean_sign, use_single_reference_value, &
4248	wall_adjustment, wall_adjustment_factor
4249
4250	IMPLICIT NONE
4251
4252	INTEGER(iwp) :: i !< loop index
4253	INTEGER(iwp) :: j !< loop index
4254	INTEGER(iwp) :: k !< loop index
4255
4256	REAL(wp) :: dvar_dz !< vertical gradient of var
4257	REAL(wp) :: l !< mixing length
4258	REAL(wp) :: l_stable !< mixing length according to stratification
4259	REAL(wp) :: ll !< adjusted l_grid
4260	REAL(wp) :: var_reference !< var at reference height
4261
4262	#if defined( __nopointer )
4263	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature
4264	#else
4265	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4266	#endif
4267
4268	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4269	IF ( dvar_dz > 0.0_wp ) THEN
4270	IF ( use_single_reference_value ) THEN
4271	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
4272	/ SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
4273	ELSE
4274	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
4275	/ SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
4276	ENDIF
4277	ELSE
4278	l_stable = l_grid(k)
4279	ENDIF
4280	!
4281	!-- Adjustment of the mixing length
4282	IF ( wall_adjustment ) THEN
4283	l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
4284	ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
4285	ELSE
4286	l = MIN( l_grid(k), l_stable )
4287	ll = l_grid(k)
4288	ENDIF
4289
4290	END SUBROUTINE mixing_length_les
4291
4292
4293	!------------------------------------------------------------------------------!
4294	! Description:
4295	! ------------
4296	!> Calculate mixing length for RANS mode.
4297	!------------------------------------------------------------------------------!
4298	SUBROUTINE mixing_length_rans( i, j, k, l, l_diss, var, var_reference )
4299
4300	USE arrays_3d, &
4301	ONLY: dd2zu
4302
4303	USE control_parameters, &
4304	ONLY: atmos_ocean_sign, use_single_reference_value
4305
4306	IMPLICIT NONE
4307
4308	INTEGER(iwp) :: i !< loop index
4309	INTEGER(iwp) :: j !< loop index
4310	INTEGER(iwp) :: k !< loop index
4311
4312	REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector
4313	REAL(wp) :: dvar_dz !< vertical gradient of var
4314	REAL(wp) :: l !< mixing length
4315	REAL(wp) :: l_diss !< mixing length for dissipation
4316	REAL(wp) :: rif !< Richardson flux number
4317	REAL(wp) :: var_reference !< var at reference height
4318
4319	#if defined( __nopointer )
4320	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature
4321	#else
4322	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4323	#endif
4324
4325	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4326
4327	duv2_dz2 = ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2 &
4328	+ ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 &
4329	+ 1E-30_wp
4330
4331	IF ( use_single_reference_value ) THEN
4332	rif = g / var_reference * dvar_dz / duv2_dz2
4333	ELSE
4334	rif = g / var(k,j,i) * dvar_dz / duv2_dz2
4335	ENDIF
4336
4337	rif = MAX( rif, -5.0_wp )
4338	rif = MIN( rif, 1.0_wp )
4339
4340	!
4341	!-- Calculate diabatic mixing length using Dyer-profile functions
4342	IF ( rif >= 0.0_wp ) THEN
4343	l = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
4344	l_diss = l
4345	ELSE
4346	!
4347	!-- In case of unstable stratification, use mixing length of neutral case
4348	!-- for l, but consider profile functions for l_diss
4349	l = l_wall(k,j,i)
4350	l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
4351	ENDIF
4352
4353	END SUBROUTINE mixing_length_rans
4354
4355
4356	!------------------------------------------------------------------------------!
4357	! Description:
4358	! ------------
4359	!> Computation of the turbulent diffusion coefficients for momentum and heat.
4360	!> @bug unstable stratification is not properly considered for kh in rans mode.
4361	!------------------------------------------------------------------------------!
4362	SUBROUTINE tcm_diffusivities( var, var_reference )
4363
4364	USE control_parameters, &
4365	ONLY: bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s, &
4366	e_min
4367
4368	USE surface_layer_fluxes_mod, &
4369	ONLY: phi_m
4370
4371	USE surface_mod, &
4372	ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, &
4373	surf_usm_h, surf_usm_v
4374
4375	INTEGER(iwp) :: i !< loop index
4376	INTEGER(iwp) :: j !< loop index
4377	INTEGER(iwp) :: k !< loop index
4378	INTEGER(iwp) :: m !< loop index
4379	INTEGER(iwp) :: n !< loop index
4380
4381	REAL(wp) :: var_reference !< reference temperature
4382
4383	#if defined( __nopointer )
4384	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature
4385	#else
4386	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4387	#endif
4388
4389	!
4390	!-- Introduce an optional minimum tke
4391	IF ( e_min > 0.0_wp ) THEN
4392	!$OMP PARALLEL DO PRIVATE(i,j,k)
4393	DO i = nxlg, nxrg
4394	DO j = nysg, nyng
4395	DO k = nzb+1, nzt
4396	e(k,j,i) = MAX( e(k,j,i), e_min ) * &
4397	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4398	ENDDO
4399	ENDDO
4400	ENDDO
4401	ENDIF
4402
4403	!
4404	!-- Call default diffusivities routine. This is always used to calculate kh.
4405	CALL tcm_diffusivities_default( var, var_reference )
4406	!
4407	!-- Call dynamic subgrid model to calculate km.
4408	IF ( les_dynamic ) THEN
4409	CALL tcm_diffusivities_dynamic
4410	ENDIF
4411
4412	!
4413	!-- In RANS mode, use MOST to calculate km and kh within the surface layer.
4414	IF ( rans_tke_e ) THEN
4415	!
4416	!-- Upward facing surfaces
4417	!-- Default surfaces
4418	n = 0
4419	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4420	DO m = 1, surf_def_h(0)%ns
4421	i = surf_def_h(0)%i(m)
4422	j = surf_def_h(0)%j(m)
4423	k = surf_def_h(0)%k(m)
4424	km(k,j,i) = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / &
4425	phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
4426	kh(k,j,i) = 1.35_wp * km(k,j,i)
4427	ENDDO
4428	!
4429	!-- Natural surfaces
4430	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4431	DO m = 1, surf_lsm_h%ns
4432	i = surf_lsm_h%i(m)
4433	j = surf_lsm_h%j(m)
4434	k = surf_lsm_h%k(m)
4435	km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / &
4436	phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
4437	kh(k,j,i) = 1.35_wp * km(k,j,i)
4438	ENDDO
4439	!
4440	!-- Urban surfaces
4441	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4442	DO m = 1, surf_usm_h%ns
4443	i = surf_usm_h%i(m)
4444	j = surf_usm_h%j(m)
4445	k = surf_usm_h%k(m)
4446	km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / &
4447	phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
4448	kh(k,j,i) = 1.35_wp * km(k,j,i)
4449	ENDDO
4450
4451	!
4452	!-- North-, south-, west and eastward facing surfaces
4453	!-- Do not consider stratification at these surfaces.
4454	DO n = 0, 3
4455	!
4456	!-- Default surfaces
4457	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4458	DO m = 1, surf_def_v(n)%ns
4459	i = surf_def_v(n)%i(m)
4460	j = surf_def_v(n)%j(m)
4461	k = surf_def_v(n)%k(m)
4462	km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m)
4463	kh(k,j,i) = 1.35_wp * km(k,j,i)
4464	ENDDO
4465	!
4466	!-- Natural surfaces
4467	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4468	DO m = 1, surf_lsm_v(n)%ns
4469	i = surf_lsm_v(n)%i(m)
4470	j = surf_lsm_v(n)%j(m)
4471	k = surf_lsm_v(n)%k(m)
4472	km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
4473	kh(k,j,i) = 1.35_wp * km(k,j,i)
4474	ENDDO
4475	!
4476	!-- Urban surfaces
4477	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4478	DO m = 1, surf_usm_v(n)%ns
4479	i = surf_usm_v(n)%i(m)
4480	j = surf_usm_v(n)%j(m)
4481	k = surf_usm_v(n)%k(m)
4482	km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
4483	kh(k,j,i) = 1.35_wp * km(k,j,i)
4484	ENDDO
4485	ENDDO
4486
4487	CALL exchange_horiz( km, nbgp )
4488	CALL exchange_horiz( kh, nbgp )
4489
4490	ENDIF
4491	!
4492	!-- Set boundary values (Neumann conditions)
4493	!-- Downward facing surfaces
4494	!$OMP PARALLEL DO PRIVATE(i,j,k)
4495	DO m = 1, bc_h(1)%ns
4496	i = bc_h(1)%i(m)
4497	j = bc_h(1)%j(m)
4498	k = bc_h(1)%k(m)
4499	km(k+1,j,i) = km(k,j,i)
4500	kh(k+1,j,i) = kh(k,j,i)
4501	ENDDO
4502	!
4503	!-- Downward facing surfaces
4504	!$OMP PARALLEL DO PRIVATE(i,j,k)
4505	DO m = 1, bc_h(0)%ns
4506	i = bc_h(0)%i(m)
4507	j = bc_h(0)%j(m)
4508	k = bc_h(0)%k(m)
4509	km(k-1,j,i) = km(k,j,i)
4510	kh(k-1,j,i) = kh(k,j,i)
4511	ENDDO
4512	!
4513	!-- Model top
4514	!$OMP PARALLEL DO
4515	DO i = nxlg, nxrg
4516	DO j = nysg, nyng
4517	km(nzt+1,j,i) = km(nzt,j,i)
4518	kh(nzt+1,j,i) = kh(nzt,j,i)
4519	ENDDO
4520	ENDDO
4521
4522	!
4523	!-- Set Neumann boundary conditions at the outflow boundaries in case of
4524	!-- non-cyclic lateral boundaries
4525	IF ( bc_radiation_l ) THEN
4526	km(:,:,nxl-1) = km(:,:,nxl)
4527	kh(:,:,nxl-1) = kh(:,:,nxl)
4528	ENDIF
4529	IF ( bc_radiation_r ) THEN
4530	km(:,:,nxr+1) = km(:,:,nxr)
4531	kh(:,:,nxr+1) = kh(:,:,nxr)
4532	ENDIF
4533	IF ( bc_radiation_s ) THEN
4534	km(:,nys-1,:) = km(:,nys,:)
4535	kh(:,nys-1,:) = kh(:,nys,:)
4536	ENDIF
4537	IF ( bc_radiation_n ) THEN
4538	km(:,nyn+1,:) = km(:,nyn,:)
4539	kh(:,nyn+1,:) = kh(:,nyn,:)
4540	ENDIF
4541
4542	END SUBROUTINE tcm_diffusivities
4543
4544
4545	!------------------------------------------------------------------------------!
4546	! Description:
4547	! ------------
4548	!> Computation of the turbulent diffusion coefficients for momentum and heat
4549	!> according to Prandtl-Kolmogorov.
4550	!------------------------------------------------------------------------------!
4551	SUBROUTINE tcm_diffusivities_default( var, var_reference )
4552
4553	USE statistics, &
4554	ONLY : rmask, sums_l_l
4555
4556	IMPLICIT NONE
4557
4558	INTEGER(iwp) :: i !< loop index
4559	INTEGER(iwp) :: j !< loop index
4560	INTEGER(iwp) :: k !< loop index
4561	!$ INTEGER(iwp) :: omp_get_thread_num !< opemmp function to get thread number
4562	INTEGER(iwp) :: sr !< statistic region
4563	INTEGER(iwp) :: tn !< thread number
4564
4565	REAL(wp) :: flag !< topography flag
4566	REAL(wp) :: l !< mixing length
4567	REAL(wp) :: ll !< adjusted mixing length
4568	REAL(wp) :: var_reference !< reference temperature
4569
4570	#if defined( __nopointer )
4571	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: var !< temperature
4572	#else
4573	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4574	#endif
4575
4576	!
4577	!-- Default thread number in case of one thread
4578	tn = 0
4579
4580	!
4581	!-- Initialization for calculation of the mixing length profile
4582	sums_l_l = 0.0_wp
4583
4584	!
4585	!-- Compute the turbulent diffusion coefficient for momentum
4586	!$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
4587	!$ tn = omp_get_thread_num()
4588
4589	IF ( les_dynamic .OR. les_mw ) THEN
4590	!$OMP DO
4591	DO i = nxlg, nxrg
4592	DO j = nysg, nyng
4593	DO k = nzb+1, nzt
4594
4595	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4596
4597	!
4598	!-- Determine the mixing length for LES closure
4599	CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
4600	!
4601	!-- Compute diffusion coefficients for momentum and heat
4602	km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
4603	kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
4604	!
4605	!-- Summation for averaged profile (cf. flow_statistics)
4606	DO sr = 0, statistic_regions
4607	sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) &
4608	* flag
4609	ENDDO
4610
4611	ENDDO
4612	ENDDO
4613	ENDDO
4614
4615	ELSEIF ( rans_tke_l ) THEN
4616
4617	!$OMP DO
4618	DO i = nxlg, nxrg
4619	DO j = nysg, nyng
4620	DO k = nzb+1, nzt
4621
4622	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4623	!
4624	!-- Mixing length for RANS mode with TKE-l closure
4625	CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
4626	!
4627	!-- Compute diffusion coefficients for momentum and heat
4628	km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
4629	kh(k,j,i) = km(k,j,i) / prandtl_number * flag
4630	!
4631	!-- Summation for averaged profile (cf. flow_statistics)
4632	DO sr = 0, statistic_regions
4633	sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) &
4634	* flag
4635	ENDDO
4636
4637	ENDDO
4638	ENDDO
4639	ENDDO
4640
4641	ELSEIF ( rans_tke_e ) THEN
4642
4643	!$OMP DO
4644	DO i = nxlg, nxrg
4645	DO j = nysg, nyng
4646	DO k = nzb+1, nzt
4647
4648	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4649	!
4650	!-- Compute diffusion coefficients for momentum and heat
4651	km(k,j,i) = c_0*4 e(k,j,i)*2 / ( diss(k,j,i) + 1.0E-30_wp ) flag
4652	kh(k,j,i) = km(k,j,i) / prandtl_number * flag
4653	!
4654	!-- Summation for averaged profile of mixing length (cf. flow_statistics)
4655	DO sr = 0, statistic_regions
4656	sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + &
4657	c_0*3 e(k,j,i) * SQRT(e(k,j,i)) / &
4658	( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
4659	ENDDO
4660
4661	ENDDO
4662	ENDDO
4663	ENDDO
4664
4665	ENDIF
4666
4667	sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundary-condition for
4668	! data output
4669	!$OMP END PARALLEL
4670
4671	END SUBROUTINE tcm_diffusivities_default
4672
4673
4674	!------------------------------------------------------------------------------!
4675	! Description:
4676	! ------------
4677	!> Calculates the eddy viscosity dynamically using the linear dynamic model
4678	!> according to
4679	!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via
4680	!> stochastic analysis."
4681	!> Monte Carlo Methods and Applications 14.4 (2008): 311-329.
4682	!>
4683	!> Furthermore dynamic bounds are used to limit the absolute value of c* as
4684	!> described in
4685	!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation:
4686	!> Stability via realizability." Physics of Fluids 29.10 (2017): 105104.
4687	!>
4688	!> @author Hauke Wurps
4689	!------------------------------------------------------------------------------!
4690	SUBROUTINE tcm_diffusivities_dynamic
4691
4692	USE arrays_3d, &
4693	ONLY: ddzw, dzw, dd2zu, w, ug, vg
4694
4695	USE grid_variables, &
4696	ONLY : ddx, ddy, dx, dy
4697
4698	IMPLICIT NONE
4699
4700	INTEGER(iwp) :: i !< running index x-direction
4701	INTEGER(iwp) :: j !< running index y-direction
4702	INTEGER(iwp) :: k !< running index z-direction
4703	INTEGER(iwp) :: l !< running index
4704	INTEGER(iwp) :: m !< running index
4705
4706	REAL(wp) :: dudx !< Gradient of u-component in x-direction
4707	REAL(wp) :: dudy !< Gradient of u-component in y-direction
4708	REAL(wp) :: dudz !< Gradient of u-component in z-direction
4709	REAL(wp) :: dvdx !< Gradient of v-component in x-direction
4710	REAL(wp) :: dvdy !< Gradient of v-component in y-direction
4711	REAL(wp) :: dvdz !< Gradient of v-component in z-direction
4712	REAL(wp) :: dwdx !< Gradient of w-component in x-direction
4713	REAL(wp) :: dwdy !< Gradient of w-component in y-direction
4714	REAL(wp) :: dwdz !< Gradient of w-component in z-direction
4715
4716	REAL(wp) :: uc(-1:1,-1:1) !< u on grid center
4717	REAL(wp) :: vc(-1:1,-1:1) !< v on grid center
4718	REAL(wp) :: wc(-1:1,-1:1) !< w on grid center
4719	! REAL(wp) :: u2(-1:1,-1:1) !< u2 on grid center
4720	! REAL(wp) :: v2(-1:1,-1:1) !< v2 on grid center
4721	! REAL(wp) :: w2(-1:1,-1:1) !< w2 on grid center
4722	! REAL(wp) :: uv(-1:1,-1:1) !< u*v on grid center
4723	! REAL(wp) :: uw(-1:1,-1:1) !< u*w on grid center
4724	! REAL(wp) :: vw(-1:1,-1:1) !< v*w on grid center
4725
4726	REAL(wp) :: ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< test filtered u
4727	REAL(wp) :: vt(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< test filtered v
4728	REAL(wp) :: wt(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< test filtered w
4729
4730	REAL(wp) :: uct !< test filtered u on grid center
4731	REAL(wp) :: vct !< test filtered v on grid center
4732	REAL(wp) :: wct !< test filtered w on grid center
4733	REAL(wp) :: u2t !< test filtered u**2 on grid center
4734	REAL(wp) :: v2t !< test filtered v**2 on grid center
4735	REAL(wp) :: w2t !< test filtered w**2 on grid center
4736	REAL(wp) :: uvt !< test filtered u*v on grid center
4737	REAL(wp) :: uwt !< test filtered u*w on grid center
4738	REAL(wp) :: vwt !< test filtered v*w on grid center
4739
4740	REAL(wp) :: sd11 !< deviatoric shear tensor
4741	REAL(wp) :: sd22 !< deviatoric shear tensor
4742	REAL(wp) :: sd33 !< deviatoric shear tensor
4743	REAL(wp) :: sd12 !< deviatoric shear tensor
4744	REAL(wp) :: sd13 !< deviatoric shear tensor
4745	REAL(wp) :: sd23 !< deviatoric shear tensor
4746
4747	REAL(wp) :: sd2 !< sum: sd_ij*sd_ij
4748
4749	REAL(wp) :: sdt11 !< filtered deviatoric shear tensor
4750	REAL(wp) :: sdt22 !< filtered deviatoric shear tensor
4751	REAL(wp) :: sdt33 !< filtered deviatoric shear tensor
4752	REAL(wp) :: sdt12 !< filtered deviatoric shear tensor
4753	REAL(wp) :: sdt13 !< filtered deviatoric shear tensor
4754	REAL(wp) :: sdt23 !< filtered deviatoric shear tensor
4755
4756	REAL(wp) :: sdt2 !< sum: sdt_ij*sdt_ij
4757
4758	REAL(wp) :: ld11 !< deviatoric stress tensor
4759	REAL(wp) :: ld22 !< deviatoric stress tensor
4760	REAL(wp) :: ld33 !< deviatoric stress tensor
4761	REAL(wp) :: ld12 !< deviatoric stress tensor
4762	REAL(wp) :: ld13 !< deviatoric stress tensor
4763	REAL(wp) :: ld23 !< deviatoric stress tensor
4764
4765	REAL(wp) :: lnn !< sum ld_nn
4766	REAL(wp) :: ldsd !< sum: ld_ij*sd_ij
4767
4768	REAL(wp) :: delta !< grid size
4769	REAL(wp) :: cst !< c*
4770	REAL(wp) :: cstnust_t !< product cnu
4771	REAL(wp) :: cst_max !< bounds of c*
4772
4773	REAL(wp), PARAMETER :: fac_cmax = 23.0_wp/(24.0_wp*sqrt(3.0_wp)) !< constant
4774
4775	!
4776	!-- velocities on grid centers:
4777	CALL tcm_box_filter_2d( u, ut )
4778	CALL tcm_box_filter_2d( v, vt )
4779	CALL tcm_box_filter_2d( w, wt )
4780
4781	DO i = nxl, nxr
4782	DO j = nys, nyn
4783	DO k = nzb+1, nzt
4784	!
4785	!-- Compute the deviatoric shear tensor s_ij on grid centers:
4786	!-- s_ij = 0.5 * ( du_i/dx_j + du_j/dx_i )
4787	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
4788	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
4789	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
4790	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
4791	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
4792
4793	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
4794	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
4795	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
4796	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
4797	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
4798
4799	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
4800	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
4801	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
4802	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
4803	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
4804
4805	sd11 = dudx
4806	sd22 = dvdy
4807	sd33 = dwdz
4808	sd12 = 0.5_wp * ( dudy + dvdx )
4809	sd13 = 0.5_wp * ( dudz + dwdx )
4810	sd23 = 0.5_wp * ( dvdz + dwdy )
4811	!
4812	!-- sum: sd_ij*sd_ij
4813	sd2 = sd112 + sd222 + sd33**2 &
4814	+ 2.0_wp * ( sd122 + sd132 + sd23**2 )
4815	!
4816	!-- The filtered velocities are needed to calculate the filtered shear
4817	!-- tensor: sdt_ij = 0.5 * ( dut_i/dx_j + dut_j/dx_i )
4818	dudx = ( ut(k,j,i+1) - ut(k,j,i) ) * ddx
4819	dudy = 0.25_wp * ( ut(k,j+1,i) + ut(k,j+1,i+1) - &
4820	ut(k,j-1,i) - ut(k,j-1,i+1) ) * ddy
4821	dudz = 0.5_wp * ( ut(k+1,j,i) + ut(k+1,j,i+1) - &
4822	ut(k-1,j,i) - ut(k-1,j,i+1) ) * dd2zu(k)
4823
4824	dvdx = 0.25_wp * ( vt(k,j,i+1) + vt(k,j+1,i+1) - &
4825	vt(k,j,i-1) - vt(k,j+1,i-1) ) * ddx
4826	dvdy = ( vt(k,j+1,i) - vt(k,j,i) ) * ddy
4827	dvdz = 0.5_wp * ( vt(k+1,j,i) + vt(k+1,j+1,i) - &
4828	vt(k-1,j,i) - vt(k-1,j+1,i) ) * dd2zu(k)
4829
4830	dwdx = 0.25_wp * ( wt(k,j,i+1) + wt(k-1,j,i+1) - &
4831	wt(k,j,i-1) - wt(k-1,j,i-1) ) * ddx
4832	dwdy = 0.25_wp * ( wt(k,j+1,i) + wt(k-1,j+1,i) - &
4833	wt(k,j-1,i) - wt(k-1,j-1,i) ) * ddy
4834	dwdz = ( wt(k,j,i) - wt(k-1,j,i) ) * ddzw(k)
4835
4836	sdt11 = dudx
4837	sdt22 = dvdy
4838	sdt33 = dwdz
4839	sdt12 = 0.5_wp * ( dudy + dvdx )
4840	sdt13 = 0.5_wp * ( dudz + dwdx )
4841	sdt23 = 0.5_wp * ( dvdz + dwdy )
4842	!
4843	!-- sum: sd_ij*sd_ij
4844	sdt2 = sdt112 + sdt222 + sdt33**2 &
4845	+ 2.0_wp * ( sdt122 + sdt132 + sdt23**2 )
4846	!
4847	!-- Need filtered velocities and filtered squared velocities on grid
4848	!-- centers. Substraction of geostrophic velocity helps to avoid
4849	!-- numerical errors in the expression <u*2> - <u><u>, which can be
4850	!-- very small (<...> means filtered).
4851	DO l = -1, 1
4852	DO m = -1, 1
4853	uc(l,m) = 0.5_wp * ( u(k,j+l,i+m) + u(k,j+l,i+m+1) ) - ug(k)
4854	vc(l,m) = 0.5_wp * ( v(k,j+l,i+m) + v(k,j+l+1,i+m) ) - vg(k)
4855	wc(l,m) = 0.5_wp * ( w(k-1,j+l,i+m) + w(k,j+l,i+m) )
4856	ENDDO
4857	ENDDO
4858
4859	CALL tcm_box_filter_2d( uc, uct )
4860	CALL tcm_box_filter_2d( vc, vct )
4861	CALL tcm_box_filter_2d( wc, wct )
4862	CALL tcm_box_filter_2d( uc**2, u2t )
4863	CALL tcm_box_filter_2d( vc**2, v2t )
4864	CALL tcm_box_filter_2d( wc**2, w2t )
4865	CALL tcm_box_filter_2d( uc*vc, uvt )
4866	CALL tcm_box_filter_2d( uc*wc, uwt )
4867	CALL tcm_box_filter_2d( vc*wc, vwt )
4868
4869	ld11 = u2t - uct*uct
4870	ld22 = v2t - vct*vct
4871	ld33 = w2t - wct*wct
4872	ld12 = uvt - uct*vct
4873	ld13 = uwt - uct*wct
4874	ld23 = vwt - vct*wct
4875
4876	lnn = ld11 + ld22 + ld33
4877	!
4878	!-- Substract trace to get deviatoric resolved stress
4879	ld11 = ld11 - lnn / 3.0_wp
4880	ld22 = ld22 - lnn / 3.0_wp
4881	ld33 = ld33 - lnn / 3.0_wp
4882
4883	ldsd = ld11 * sdt11 + ld22 * sdt22 + ld33 * sdt33 + &
4884	2.0_wp * ( ld12 * sdt12 + ld13 * sdt13 + ld23 * sdt23 )
4885	!
4886	!-- c* nu*^T is SGS viscosity on test filter level:
4887	cstnust_t = -ldsd / sdt2
4888	!
4889	!-- The model was only tested for an isotropic grid. The following
4890	!-- expression was a recommendation of Stefan Heinz.
4891	delta = MAX( dx, dy, dzw(k) )
4892	cst = cstnust_t / ( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) )
4893	!
4894	!-- Calculate border according to Mokhtarpoor and Heinz (2017)
4895	cst_max = fac_cmax * SQRT( e(k,j,i) ) / ( delta * SQRT( 2.0_wp * sd2 ) )
4896
4897	IF ( ABS( cst ) > cst_max ) THEN
4898	cst = cst_max * cst / ABS( cst )
4899	ENDIF
4900
4901	km(k,j,i) = cst * delta * SQRT( e(k,j,i) )
4902
4903	ENDDO
4904	ENDDO
4905	ENDDO
4906
4907	END SUBROUTINE tcm_diffusivities_dynamic
4908
4909
4910	!------------------------------------------------------------------------------!
4911	! Description:
4912	! ------------
4913	!> This subroutine acts as a box filter with filter width 2 * dx.
4914	!> Output is only one point.
4915	!------------------------------------------------------------------------------!
4916	SUBROUTINE tcm_box_filter_2d_single( var, var_fil )
4917
4918	IMPLICIT NONE
4919
4920	REAL(wp) :: var(-1:1,-1:1) !< variable to be filtered
4921	REAL(wp) :: var_fil !< filtered variable
4922	!
4923	!-- It is assumed that a box with a side length of 2 * dx and centered at the
4924	!-- variable*s position contains one half of the four closest neigbours and one
4925	!-- forth of the diagonally closest neighbours.
4926	var_fil = 0.25_wp * ( var(0,0) + &
4927	0.5_wp * ( var(0,1) + var(1,0) + &
4928	var(0,-1) + var(-1,0) ) + &
4929	0.25_wp * ( var(1,1) + var(1,-1) + &
4930	var(-1,1) + var(-1,-1) ) )
4931
4932	END SUBROUTINE tcm_box_filter_2d_single
4933
4934	!------------------------------------------------------------------------------!
4935	! Description:
4936	! ------------
4937	!> This subroutine acts as a box filter with filter width 2 * dx.
4938	!> The filtered variable var_fil is on the same grid as var.
4939	!------------------------------------------------------------------------------!
4940	SUBROUTINE tcm_box_filter_2d_array( var, var_fil )
4941
4942	IMPLICIT NONE
4943
4944	INTEGER(iwp) :: i !< running index x-direction
4945	INTEGER(iwp) :: j !< running index y-direction
4946	INTEGER(iwp) :: k !< running index z-direction
4947
4948	REAL(wp) :: var(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< variable to be filtered
4949	REAL(wp) :: var_fil(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< filtered variable
4950	!
4951	!-- It is assumed that a box with a side length of 2 * dx and centered at the
4952	!-- variable's position contains one half of the four closest neigbours and one
4953	!-- forth of the diagonally closest neighbours.
4954	DO i = nxlg+1, nxrg-1
4955	DO j = nysg+1, nyng-1
4956	DO k = nzb, nzt+1
4957	var_fil(k,j,i) = 0.25_wp * ( var(k,j,i) + &
4958	0.5_wp * ( var(k,j,i+1) + var(k,j+1,i) + &
4959	var(k,j,i-1) + var(k,j-1,i) ) +&
4960	0.25_wp * ( var(k,j+1,i+1) + var(k,j+1,i-1) + &
4961	var(k,j-1,i+1) + var(k,j-1,i-1) ) )
4962	END DO
4963	END DO
4964	END DO
4965
4966	END SUBROUTINE tcm_box_filter_2d_array
4967
4968
4969	!------------------------------------------------------------------------------!
4970	! Description:
4971	! ------------
4972	!> Swapping of timelevels.
4973	!------------------------------------------------------------------------------!
4974	SUBROUTINE tcm_swap_timelevel ( mod_count )
4975
4976	IMPLICIT NONE
4977
4978
4979	#if defined( __nopointer )
4980	INTEGER(iwp) :: i !< loop index x direction
4981	INTEGER(iwp) :: j !< loop index y direction
4982	INTEGER(iwp) :: k !< loop index z direction
4983	#endif
4984	INTEGER, INTENT(IN) :: mod_count !< flag defining where pointers point to
4985
4986	#if defined( __nopointer )
4987
4988	IF ( .NOT. constant_diffusion ) THEN
4989	DO i = nxlg, nxrg
4990	DO j = nysg, nyng
4991	DO k = nzb, nzt+1
4992	e(k,j,i) = e_p(k,j,i)
4993	ENDDO
4994	ENDDO
4995	ENDDO
4996	ENDIF
4997
4998	IF ( rans_tke_e ) THEN
4999	DO i = nxlg, nxrg
5000	DO j = nysg, nyng
5001	DO k = nzb, nzt+1
5002	diss(k,j,i) = diss_p(k,j,i)
5003	ENDDO
5004	ENDDO
5005	ENDDO
5006	ENDIF
5007
5008	#else
5009
5010	SELECT CASE ( mod_count )
5011
5012	CASE ( 0 )
5013
5014	IF ( .NOT. constant_diffusion ) THEN
5015	e => e_1; e_p => e_2
5016	ENDIF
5017
5018	IF ( rans_tke_e ) THEN
5019	diss => diss_1; diss_p => diss_2
5020	ENDIF
5021
5022	CASE ( 1 )
5023
5024	IF ( .NOT. constant_diffusion ) THEN
5025	e => e_2; e_p => e_1
5026	ENDIF
5027
5028	IF ( rans_tke_e ) THEN
5029	diss => diss_2; diss_p => diss_1
5030	ENDIF
5031
5032	END SELECT
5033	#endif
5034
5035	END SUBROUTINE tcm_swap_timelevel
5036
5037
5038	END MODULE turbulence_closure_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |