Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |