Home

Context Navigation

turbulence_closure_mod.f90 @ 3014

Last change on this file since 3014 was 3014, checked in by maronga, 5 years ago
series of bugfixes
Property svn:keywords set to `Id`
File size: 183.6 KB

Line
1	!> @file turbulence_closure_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 2017-2018 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: turbulence_closure_mod.f90 3014 2018-05-09 08:42:38Z maronga $
27	! Bugfix: nzb_do and nzt_do were not used for 3d data output
28	!
29	! 3004 2018-04-27 12:33:25Z Giersch
30	! Further allocation checks implemented
31	!
32	! 2938 2018-03-27 15:52:42Z suehring
33	! Further todo's
34	!
35	! 2936 2018-03-27 14:49:27Z suehring
36	! - defined l_grid only within this module
37	! - Moved l_wall definition from modules.f90
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
43	! - Moved init_mixing_length from init_grid.f90 and
44	! renamed it to tcm_init_mixing_length
45	!
46	! 2764 2018-01-22 09:25:36Z gronemeier
47	! Bugfix: remove duplicate SAVE statements
48	!
49	! 2746 2018-01-15 12:06:04Z suehring
50	! Move flag plant canopy to modules
51	!
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
60	!
61	! 2696 2017-12-14 17:12:51Z kanani
62	! Initial revision
63	!
64	!
65	!
66	!
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
78	!> add OpenMP directives whereever possible
79	!> remove debug output variables (dummy1, dummy2, dummy3)
80	!> @todo Move initialization of wall-mixing length from init_grid
81	!> @todo Check for random disturbances
82	!> @note <Enter notes on the module>
83	!> @bug TKE-e closure still crashes due to too small dt
84	!------------------------------------------------------------------------------!
85	MODULE turbulence_closure_mod
86
87
88	#if defined( __nopointer )
89	USE arrays_3d, &
90	ONLY: diss, diss_p, dzu, e, e_p, kh, km, &
91	mean_inflow_profiles, prho, pt, tdiss_m, te_m, tend, u, v, vpt, w
92	#else
93	USE arrays_3d, &
94	ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, &
95	e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, &
96	te_m, tend, u, v, vpt, w
97	#endif
98
99	USE control_parameters, &
100	ONLY: constant_diffusion, dt_3d, e_init, humidity, inflow_l, &
101	initializing_actions, intermediate_timestep_count, &
102	intermediate_timestep_count_max, kappa, km_constant, les_mw, &
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
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, &
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
127
128	USE kinds
129
130	USE pegrid
131
132	USE plant_canopy_model_mod, &
133	ONLY: pcm_tendency
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
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
153	REAL(wp) :: sig_diss = 1.3_wp !< factor to calculate K_diss from Km
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
156
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
159
160	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< near-wall mixing length
161
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
165
166
167	PUBLIC c_m, dummy1, dummy2, dummy3
168
169	!
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
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
181
182	!
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
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
211
212	!
213	!-- Initialization of arrays
214	INTERFACE tcm_init_arrays
215	MODULE PROCEDURE tcm_init_arrays
216	END INTERFACE tcm_init_arrays
217
218	!
219	!-- Initialization of TKE production term
220	INTERFACE production_e_init
221	MODULE PROCEDURE production_e_init
222	END INTERFACE production_e_init
223
224	!
225	!-- Prognostic equations for TKE and TKE dissipation rate
226	INTERFACE tcm_prognostic
227	MODULE PROCEDURE tcm_prognostic
228	MODULE PROCEDURE tcm_prognostic_ij
229	END INTERFACE tcm_prognostic
230
231	!
232	!-- Production term for TKE
233	INTERFACE production_e
234	MODULE PROCEDURE production_e
235	MODULE PROCEDURE production_e_ij
236	END INTERFACE production_e
237
238	!
239	!-- Diffusion term for TKE
240	INTERFACE diffusion_e
241	MODULE PROCEDURE diffusion_e
242	MODULE PROCEDURE diffusion_e_ij
243	END INTERFACE diffusion_e
244
245	!
246	!-- Diffusion term for TKE dissipation rate
247	INTERFACE diffusion_diss
248	MODULE PROCEDURE diffusion_diss
249	MODULE PROCEDURE diffusion_diss_ij
250	END INTERFACE diffusion_diss
251
252	!
253	!-- Mixing length for LES case
254	INTERFACE mixing_length_les
255	MODULE PROCEDURE mixing_length_les
256	END INTERFACE mixing_length_les
257
258	!
259	!-- Mixing length for RANS case
260	INTERFACE mixing_length_rans
261	MODULE PROCEDURE mixing_length_rans
262	END INTERFACE mixing_length_rans
263
264	!
265	!-- Calculate diffusivities
266	INTERFACE tcm_diffusivities
267	MODULE PROCEDURE tcm_diffusivities
268	END INTERFACE tcm_diffusivities
269
270	!
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
275
276	SAVE
277
278	PRIVATE
279	!
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
285
286
287	CONTAINS
288
289	!------------------------------------------------------------------------------!
290	! Description:
291	! ------------
292	!> Check parameters routine for turbulence closure module.
293	!------------------------------------------------------------------------------!
294	SUBROUTINE tcm_check_parameters
295
296	USE control_parameters, &
297	ONLY: message_string, nest_domain, neutral, turbulent_inflow, &
298	turbulent_outflow
299
300	IMPLICIT NONE
301
302	!
303	!-- Define which turbulence closure is going to be used
304	IF ( rans_mode ) THEN
305
306	c_m = 0.4_wp !according to Detering and Etling (1985)
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
316	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) == 0 &
317	.AND. .NOT. nest_domain ) THEN
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
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
354	IF ( rans_tke_e ) THEN
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	!------------------------------------------------------------------------------!
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' )
396	unit = 'mixing length'
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
428	!
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	!------------------------------------------------------------------------------!
464	! Description:
465	! ------------
466	!> Average 3D data.
467	!------------------------------------------------------------------------------!
468	SUBROUTINE tcm_3d_data_averaging( mode, variable )
469
470
471	USE averaging, &
472	ONLY: diss_av, kh_av, km_av
473
474	USE control_parameters, &
475	ONLY: average_count_3d
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
492	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
493	ENDIF
494	diss_av = 0.0_wp
495
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
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' )
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
524	ENDDO
525	ENDDO
526	ENDIF
527
528	CASE ( 'kh' )
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
535	ENDDO
536	ENDDO
537	ENDIF
538
539	CASE ( 'km' )
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
546	ENDDO
547	ENDDO
548	ENDIF
549
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' )
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
567	ENDDO
568	ENDDO
569	ENDIF
570
571	CASE ( 'kh' )
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
579	ENDDO
580	ENDDO
581	ENDIF
582
583	CASE ( 'km' )
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
591	ENDDO
592	ENDDO
593	ENDIF
594
595	END SELECT
596
597	ENDIF
598
599	END SUBROUTINE tcm_3d_data_averaging
600
601
602	!------------------------------------------------------------------------------!
603	! Description:
604	! ------------
605	!> Define 2D output variables.
606	!------------------------------------------------------------------------------!
607	SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, &
608	two_d, nzb_do, nzt_do )
609
610	USE averaging, &
611	ONLY: diss_av, kh_av, km_av
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
629	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
630
631	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
632	!< array to which output data is resorted to
633
634	found = .TRUE.
635
636	SELECT CASE ( TRIM( variable ) )
637
638
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
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
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
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
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
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
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
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	! ------------
726	!> Define 3D output variables.
727	!------------------------------------------------------------------------------!
728	SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
729
730
731	USE averaging, &
732	ONLY: diss_av, kh_av, km_av
733
734	IMPLICIT NONE
735
736	CHARACTER (LEN=*) :: variable !<
737
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)
744
745	LOGICAL :: found !<
746
747	REAL(wp) :: fill_value = -999.0_wp !< value for the _FillValue attribute
748
749	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
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
761	DO i = nxl, nxr
762	DO j = nys, nyn
763	DO k = nzb_do, nzt_do
764	local_pf(i,j,k) = diss(k,j,i)
765	ENDDO
766	ENDDO
767	ENDDO
768	ELSE
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
773	DO i = nxl, nxr
774	DO j = nys, nyn
775	DO k = nzb_do, nzt_do
776	local_pf(i,j,k) = diss_av(k,j,i)
777	ENDDO
778	ENDDO
779	ENDDO
780	ENDIF
781
782	CASE ( 'kh' )
783	IF ( av == 0 ) THEN
784	DO i = nxl, nxr
785	DO j = nys, nyn
786	DO k = nzb_do, nzt_do
787	local_pf(i,j,k) = kh(k,j,i)
788	ENDDO
789	ENDDO
790	ENDDO
791	ELSE
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
796	DO i = nxl, nxr
797	DO j = nys, nyn
798	DO k = nzb_do, nzt_do
799	local_pf(i,j,k) = kh_av(k,j,i)
800	ENDDO
801	ENDDO
802	ENDDO
803	ENDIF
804
805	CASE ( 'km' )
806	IF ( av == 0 ) THEN
807	DO i = nxl, nxr
808	DO j = nys, nyn
809	DO k = nzb_do, nzt_do
810	local_pf(i,j,k) = km(k,j,i)
811	ENDDO
812	ENDDO
813	ENDDO
814	ELSE
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
819	DO i = nxl, nxr
820	DO j = nys, nyn
821	DO k = nzb_do, nzt_do
822	local_pf(i,j,k) = km_av(k,j,i)
823	ENDDO
824	ENDDO
825	ENDDO
826	ENDIF
827
828	CASE ( 'dummy1' ) !### remove later
829	IF ( av == 0 ) THEN
830	DO i = nxl, nxr
831	DO j = nys, nyn
832	DO k = nzb_do, nzt_do
833	local_pf(i,j,k) = dummy1(k,j,i)
834	ENDDO
835	ENDDO
836	ENDDO
837	ENDIF
838
839	CASE ( 'dummy2' ) !### remove later
840	IF ( av == 0 ) THEN
841	DO i = nxl, nxr
842	DO j = nys, nyn
843	DO k = nzb_do, nzt_do
844	local_pf(i,j,k) = dummy2(k,j,i)
845	ENDDO
846	ENDDO
847	ENDDO
848	ENDIF
849
850	CASE ( 'dummy3' ) !### remove later
851	IF ( av == 0 ) THEN
852	DO i = nxl, nxr
853	DO j = nys, nyn
854	DO k = nzb_do, nzt_do
855	local_pf(i,j,k) = dummy3(k,j,i)
856	ENDDO
857	ENDDO
858	ENDDO
859	ENDIF
860
861	CASE DEFAULT
862	found = .FALSE.
863
864	END SELECT
865
866	END SUBROUTINE tcm_data_output_3d
867
868
869	!------------------------------------------------------------------------------!
870	! Description:
871	! ------------
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
882	USE pmc_interface, &
883	ONLY: nested_run
884
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
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.
911	IF ( rans_tke_e .OR. use_sgs_for_particles .OR. wang_kernel .OR. &
912	collision_turbulence .OR. nested_run ) THEN
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) )
921	IF ( rans_tke_e .OR. nested_run ) THEN
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. &
934	wang_kernel .OR. collision_turbulence .OR. nested_run ) THEN
935	diss => diss_1
936	IF ( rans_tke_e .OR. nested_run ) THEN
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	! ------------
948	!> Initialization of turbulence closure module.
949	!------------------------------------------------------------------------------!
950	SUBROUTINE tcm_init
951
952	USE control_parameters, &
953	ONLY: complex_terrain, dissipation_1d, topography
954
955	USE model_1d_mod, &
956	ONLY: diss1d, e1d, kh1d, km1d, l1d
957
958	USE surface_mod, &
959	ONLY: get_topography_top_index_ji
960
961	IMPLICIT NONE
962
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 !<
968
969	!
970	!-- Initialize mixing length
971	CALL tcm_init_mixing_length
972
973	!
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
990	e = 0.0_wp
991	ENDIF
992
993	IF ( rans_tke_e ) THEN
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?
997	diss(:,j,i) = diss1d
998	ENDDO
999	ENDDO
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
1009	ENDIF
1010
1011	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 .OR. &
1012	INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN
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
1030	! production terms, as long as not yet
1031	! e = (u/cm)*2 at k=nzb+1
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
1048	IF ( rans_tke_e ) THEN
1049	tdiss_m = 0.0_wp
1050	diss_p = diss
1051	ENDIF
1052
1053	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &
1054	TRIM( initializing_actions ) == 'cyclic_fill' ) &
1055	THEN
1056
1057	!
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	!
1075	!-- Initialization of the turbulence recycling method
1076	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
1077	turbulent_inflow ) THEN
1078	mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e
1079	!
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	!
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	!
1112	!-- Inside buildings set TKE back to zero.
1113	!-- Other scalars (km, kh, diss, ...) are ignored at present,
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
1142	! Description:
1143	! -----------------------------------------------------------------------------!
1144	!> Pre-computation of grid-dependent and near-wall mixing length.
1145	!------------------------------------------------------------------------------!
1146	SUBROUTINE tcm_init_mixing_length
1147
1148	USE arrays_3d, &
1149	ONLY: dzw, ug, vg, zu, zw
1150
1151	USE control_parameters, &
1152	ONLY: bc_lr_cyc, bc_ns_cyc, f, kappa, message_string, &
1153	wall_adjustment_factor
1154
1155	USE grid_variables, &
1156	ONLY: dx, dy
1157
1158	USE indices, &
1159	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, &
1160	nzt, wall_flags_0
1161
1162	USE kinds
1163
1164
1165	IMPLICIT NONE
1166
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
1184
1185	INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
1186
1187	INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
1188
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
1192	REAL(wp) :: radius !< search radius in meter
1193
1194	ALLOCATE( l_grid(1:nzt) )
1195	ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1196	!
1197	!-- Initialize the mixing length in case of an LES-simulation
1198	IF ( .NOT. rans_mode ) THEN
1199	!
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)
1212
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
1225	!
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
1233	!
1234	!-- Check if current gridpoint belongs to the atmosphere
1235	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1236	!
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) )
1241	!
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) )
1245	!
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 )
1250	!
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 )
1255	!
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 ) )
1262	!
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 ) )
1269	!
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 ) )
1276	!
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 ) )
1283	!
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 ) ) )
1291	!
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 ) )
1300	!
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
1312	ENDDO
1313	ENDDO
1314
1315	ELSE
1316	!
1317	!-- Initialize the mixing length in case of a RANS simulation
1318
1319	!
1320	!-- Calculate mixing length according to Blackadar (1962)
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
1334	!
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	!
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	!
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
1378
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
1393	IF ( zu(k+kk)-zu(k) >= radius .OR. k+kk >= k_max_topo ) EXIT
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
1408	IF ( zu(k-rad_k_b) > zu(k_max_topo) ) rad_k_b = 0
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) )
1420	ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
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	!
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 )
1438
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 )
1446	!
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
1499
1500	IF ( bc_ns_cyc ) THEN
1501	!
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	!
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 (+++)
1529	vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1530	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1531	shortest_distance( vic_yz, .TRUE., ii ) )
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").
1537	vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1538	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1539	shortest_distance( vic_yz, .TRUE., ii ) )
1540
1541	ENDIF
1542	!
1543	!-- Search for walls within octant (+--)
1544	vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1545	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1546	shortest_distance( vic_yz, .FALSE., ii ) )
1547	!
1548	!-- Search for walls within octant (++-)
1549	vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1550	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1551	shortest_distance( vic_yz, .FALSE., ii ) )
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 (-++)
1566	vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1567	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1568	shortest_distance( vic_yz, .TRUE., -ii ) )
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").
1574	vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1575	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1576	shortest_distance( vic_yz, .TRUE., -ii ) )
1577
1578	ENDIF
1579	!
1580	!-- Search for walls within octant (---)
1581	vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1582	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1583	shortest_distance( vic_yz, .FALSE., -ii ) )
1584	!
1585	!-- Search for walls within octant (-+-)
1586	vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1587	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1588	shortest_distance( vic_yz, .FALSE., -ii ) )
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
1606	DEALLOCATE( vicinity )
1607	DEALLOCATE( vic_yz )
1608
1609	ENDIF !check vertical size of vicinity
1610
1611	ENDDO !k loop
1612
1613	DEALLOCATE( wall_flags_0_global )
1614
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
1639	INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
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
1648
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
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.
1669	DO jj = 0, rad_j
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	) &
1677	)
1678	ENDDO
1679	ENDIF
1680
1681	END FUNCTION
1682
1683	!------------------------------------------------------------------------------!
1684	! Description:
1685	! ------------
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.
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
1709	DO i = il, ir
1710	DO j = js, jn
1711	DO k = kb, kt
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
1720	END SUBROUTINE tcm_init_mixing_length
1721
1722
1723	!------------------------------------------------------------------------------!
1724	! Description:
1725	! ------------
1726	!> Initialize virtual velocities used later in production_e.
1727	!------------------------------------------------------------------------------!
1728	SUBROUTINE production_e_init
1729
1730	USE arrays_3d, &
1731	ONLY: drho_air_zw, zu
1732
1733	USE control_parameters, &
1734	ONLY: constant_flux_layer
1735
1736	USE surface_mod, &
1737	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h
1738
1739	IMPLICIT NONE
1740
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
1745
1746	IF ( constant_flux_layer ) THEN
1747	!
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
1761
1762	i = surf_def_h(0)%i(m)
1763	j = surf_def_h(0)%j(m)
1764	k = surf_def_h(0)%k(m)
1765	!
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 )
1778
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)
1782
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
1788	!
1789	!-- Default surfaces, downward-facing surfaces
1790	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1791	DO m = 1, surf_def_h(1)%ns
1792
1793	i = surf_def_h(1)%i(m)
1794	j = surf_def_h(1)%j(m)
1795	k = surf_def_h(1)%k(m)
1796
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 )
1805
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)
1809
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)
1813
1814	ENDDO
1815	!
1816	!-- Natural surfaces, upward-facing
1817	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1818	DO m = 1, surf_lsm_h%ns
1819
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 )
1836
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
1846	!
1847	!-- Urban surfaces, upward-facing
1848	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1849	DO m = 1, surf_usm_h%ns
1850
1851	i = surf_usm_h%i(m)
1852	j = surf_usm_h%j(m)
1853	k = surf_usm_h%k(m)
1854	!
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 )
1867
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)
1871
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)
1875
1876	ENDDO
1877
1878	ENDIF
1879
1880	END SUBROUTINE production_e_init
1881
1882
1883	!------------------------------------------------------------------------------!
1884	! Description:
1885	! ------------
1886	!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
1887	!> Vector-optimized version
1888	!------------------------------------------------------------------------------!
1889	SUBROUTINE tcm_prognostic
1890
1891	USE arrays_3d, &
1892	ONLY: ddzu
1893
1894	USE control_parameters, &
1895	ONLY: f, scalar_advec, tsc
1896
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
1900
1901	IMPLICIT NONE
1902
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
1909
1910	REAL(wp) :: sbt !< wheighting factor for sub-time step
1911
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
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
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
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	!
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	!
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
2041	ENDIF ! TKE equation
2042
2043	!
2044	!-- If required, compute prognostic equation for TKE dissipation rate
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
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
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	!
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	!
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	! ------------
2172	!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2173	!> Cache-optimized version
2174	!------------------------------------------------------------------------------!
2175	SUBROUTINE tcm_prognostic_ij( i, j, i_omp, tn )
2176
2177	USE arrays_3d, &
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
2180
2181	USE control_parameters, &
2182	ONLY: f, tsc
2183
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
2188	IMPLICIT NONE
2189
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 !<
2198
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
2201
2202	!
2203	!-- If required, compute prognostic equation for turbulent kinetic
2204	!-- energy (TKE)
2205	IF ( .NOT. constant_diffusion ) THEN
2206
2207	!
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
2221
2222	advec(:) = tend(:,j,i)
2223
2224	CALL production_e( i, j )
2225
2226	produc(:) = tend(:,j,i) - advec(:)
2227
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
2237
2238	!
2239	!-- Additional sink term for flows through plant canopies
2240	IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 )
2241
2242	CALL user_actions( i, j, 'e-tendency' )
2243
2244	!
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
2258
2259	!
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
2269
2270	!
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
2285
2286	ENDIF ! TKE equation
2287
2288	!
2289	!-- If required, compute prognostic equation for TKE dissipation rate
2290	IF ( rans_tke_e ) THEN
2291
2292	!
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
2306
2307	!
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
2315
2316	CALL diffusion_diss( i, j )
2317
2318	!
2319	!-- Additional sink term for flows through plant canopies
2320	! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !### not yet implemented
2321
2322	! CALL user_actions( i, j, 'diss-tendency' ) !### not yet implemented
2323
2324	!
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) ) &
2332	) &
2333	* MERGE( 1.0_wp, 0.0_wp, &
2334	BTEST( wall_flags_0(k,j,i), 0 )&
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
2338
2339	!
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
2347
2348	!
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
2363
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)
2367
2368	ENDIF ! dissipation equation
2369
2370	END SUBROUTINE tcm_prognostic_ij
2371
2372
2373	!------------------------------------------------------------------------------!
2374	! Description:
2375	! ------------
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!
2380	!------------------------------------------------------------------------------!
2381	SUBROUTINE production_e
2382
2383	USE arrays_3d, &
2384	ONLY: ddzw, dd2zu, drho_air_zw, q, ql
2385
2386	USE cloud_parameters, &
2387	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
2388
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
2393
2394	USE grid_variables, &
2395	ONLY: ddx, dx, ddy, dy
2396
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
2400
2401	IMPLICIT NONE
2402
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
2410
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"
2423
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
2433
2434	DO i = nxl, nxr
2435
2436	IF ( constant_flux_layer ) THEN
2437
2438	!
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
2444
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)
2451
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)
2458
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)
2464
2465	ENDDO
2466	ENDDO
2467
2468	!
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
2473	!
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
2487	!
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
2494	!
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
2505	!
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
2512	!
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
2523	!
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
2531	!
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)
2540
2541	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2542	* 0.5_wp * dx
2543	!
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
2550	!
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)
2558
2559	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2560	* 0.5_wp * dx
2561	!
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
2568	!
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)
2576
2577	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2578	* 0.5_wp * dx
2579	!
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
2587	!
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)
2593	!
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
2603	!
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)
2609
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
2614	!
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)
2620
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
2625	!
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)
2632
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)
2635
2636	ENDDO
2637	ENDDO
2638
2639	DO j = nys, nyn
2640	DO k = nzb+1, nzt
2641
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) )
2647
2648	IF ( def < 0.0_wp ) def = 0.0_wp
2649
2650	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
2651
2652	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
2653
2654	ENDDO
2655	ENDDO
2656
2657	ELSE
2658
2659	DO j = nys, nyn
2660	!
2661	!-- Calculate TKE production by shear. Here, no additional
2662	!-- wall-bounded code is considered.
2663	!-- Why?
2664	DO k = nzb+1, nzt
2665
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)
2672
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)
2679
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	)
2696
2697	IF ( def < 0.0_wp ) def = 0.0_wp
2698
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
2703	ENDDO
2704	ENDDO
2705
2706	ENDIF
2707
2708	!
2709	!-- If required, calculate TKE production by buoyancy
2710	IF ( .NOT. neutral ) THEN
2711
2712	IF ( .NOT. humidity ) THEN
2713
2714	IF ( ocean ) THEN
2715	!
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
2733	!
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)
2740	DO m = surf_s, surf_e
2741	k = surf_def_h(l)%k(m)
2742	tend(k,j,i) = tend(k,j,i) + g / &
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)
2747	ENDDO
2748	ENDDO
2749
2750	ENDIF
2751
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)
2762	ENDDO
2763	ENDIF
2764
2765	ENDDO
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
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	)
2786	ENDDO
2787
2788	IF ( use_surface_fluxes ) THEN
2789	!
2790	!-- Default surfaces, up- and downward-facing
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)
2796	tend(k,j,i) = tend(k,j,i) + g / &
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
2802	ENDDO
2803	!
2804	!-- Natural surfaces
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)
2809	tend(k,j,i) = tend(k,j,i) + g / &
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)
2814	ENDDO
2815	!
2816	!-- Urban surfaces
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
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
2828
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)
2834	tend(k,j,i) = tend(k,j,i) + g / &
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)
2839	ENDDO
2840	ENDIF
2841	ENDDO
2842
2843	ENDIF
2844
2845	ELSE
2846
2847	DO j = nys, nyn
2848
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
2914
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)
2927	DO m = surf_s, surf_e
2928	k = surf_def_h(l)%k(m)
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
2952	tend(k,j,i) = tend(k,j,i) + g / &
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)
2958	ENDDO
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)
2966
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)
2995	ENDDO
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)
3002
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
3024
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
3035	ENDIF
3036
3037	IF ( use_top_fluxes ) THEN
3038
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
3077	ENDDO
3078
3079	ENDIF
3080
3081	ENDIF
3082
3083	ENDIF
3084
3085	ENDDO
3086
3087	END SUBROUTINE production_e
3088
3089
3090	!------------------------------------------------------------------------------!
3091	! Description:
3092	! ------------
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!
3097	!------------------------------------------------------------------------------!
3098	SUBROUTINE production_e_ij( i, j )
3099
3100	USE arrays_3d, &
3101	ONLY: ddzw, dd2zu, drho_air_zw, q, ql
3102
3103	USE cloud_parameters, &
3104	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
3105
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
3110
3111	USE grid_variables, &
3112	ONLY: ddx, dx, ddy, dy
3113
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
3117
3118	IMPLICIT NONE
3119
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
3127
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"
3140
3141
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
3151
3152
3153	IF ( constant_flux_layer ) THEN
3154	!
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
3159
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)
3165
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)
3171
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)
3177
3178	ENDDO
3179	!
3180	!-- Compute gradients at north- and south-facing surfaces.
3181	!-- Note, no vertical natural surfaces so far.
3182	DO l = 0, 1
3183	!
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)
3191
3192	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3193	* 0.5_wp * dy
3194	!
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
3201	!
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)
3209
3210	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3211	* 0.5_wp * dy
3212	!
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
3219	!
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)
3227
3228	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3229	* 0.5_wp * dy
3230	!
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
3238	!
3239	!-- Compute gradients at east- and west-facing walls
3240	DO l = 2, 3
3241	!
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)
3249
3250	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3251	* 0.5_wp * dx
3252	!
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
3259	!
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)
3267
3268	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3269	* 0.5_wp * dx
3270	!
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
3277	!
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)
3285
3286	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3287	* 0.5_wp * dx
3288	!
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
3296	!
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)
3303	!
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)
3311
3312	ENDDO
3313	!
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)
3319
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
3323	!
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)
3329
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
3333	!
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)
3340
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)
3343
3344	ENDDO
3345
3346	DO k = nzb+1, nzt
3347
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) )
3352
3353	!
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	! )
3365
3366	IF ( def < 0.0_wp ) def = 0.0_wp
3367
3368	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3369
3370	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
3371
3372	ENDDO
3373
3374	ELSE
3375	!
3376	!-- Calculate TKE production by shear. Here, no additional
3377	!-- wall-bounded code is considered.
3378	!-- Why?
3379	DO k = nzb+1, nzt
3380
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)
3386
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)
3392
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)
3398
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) )
3403
3404	!
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	! )
3416
3417	IF ( def < 0.0_wp ) def = 0.0_wp
3418
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
3422
3423	ENDDO
3424
3425	ENDIF
3426
3427	!
3428	!-- If required, calculate TKE production by buoyancy
3429	IF ( .NOT. neutral ) THEN
3430
3431	IF ( .NOT. humidity ) THEN
3432
3433	IF ( ocean ) THEN
3434	!
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
3438
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
3452
3453	IF ( use_surface_fluxes ) THEN
3454	!
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)
3459	DO m = surf_s, surf_e
3460	k = surf_def_h(l)%k(m)
3461	tend(k,j,i) = tend(k,j,i) + g / &
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)
3466	ENDDO
3467	ENDDO
3468
3469	ENDIF
3470
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)
3481	ENDDO
3482	ENDIF
3483
3484	ELSE
3485
3486	DO k = nzb+1, nzt
3487	!
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	)
3501
3502	ENDDO
3503
3504	IF ( use_surface_fluxes ) THEN
3505	!
3506	!-- Default surfaces, up- and downward-facing
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)
3512	tend(k,j,i) = tend(k,j,i) + g / &
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)
3517	ENDDO
3518	ENDDO
3519	!
3520	!-- Natural surfaces
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)
3525	tend(k,j,i) = tend(k,j,i) + g / &
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)
3530	ENDDO
3531	!
3532	!-- Urban surfaces
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)
3537	tend(k,j,i) = tend(k,j,i) + g / &
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)
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)
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
3557
3558	ENDIF
3559
3560	ELSE
3561
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
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
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 )
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
3657	tend(k,j,i) = tend(k,j,i) + g / &
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)
3663	ENDDO
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)
3671
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
3693
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
3738	ENDIF
3739
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