Home

Context Navigation

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

Last change on this file since 3602 was 3550, checked in by gronemeier, 5 years ago
calculate diss production same in vector and cache optimization; move boundary condition for e and diss to boundary_conds
Property svn:keywords set to `Id`
File size: 186.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |