Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |