Home

Context Navigation

turbulence_closure_mod.f90 @ 3300

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