Home

Context Navigation

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

Last change on this file since 3386 was 3386, checked in by gronemeier, 7 years ago
renamed tcm_prognostic to tcm_prognostic_equations
Property svn:keywords set to `Id`
File size: 209.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |