Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |