Home

Context Navigation

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

Last change on this file since 2809 was 2764, checked in by gronemeier, 7 years ago
Bugfix: remove duplicate SAVE statements
Property svn:keywords set to `Id`
File size: 153.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |