Home

Context Navigation

turbulence_closure_mod.f90 @ 3359

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