Home

Context Navigation

source: palm/tags/release-5.0/SOURCE/turbulence_closure_mod.f90 @ 4106

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |