Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |