Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |