Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |