Home

Context Navigation

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

Last change on this file since 3646 was 3646, checked in by kanani, 6 years ago
Bugfix: replace simulated_time by time_since_reference_point where required
Property svn:keywords set to `Id`
File size: 193.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-2018 Leibniz Universitaet Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: turbulence_closure_mod.f90 3646 2018-12-28 17:58:49Z kanani $
27	! Remove unused variable simulated_time
28	!
29	! 3636 2018-12-19 13:48:34Z raasch
30	! nopointer option removed
31	!
32	! 3634 2018-12-18 12:31:28Z knoop
33	! OpenACC port for SPEC
34	!
35	! 3550 2018-11-21 16:01:01Z gronemeier
36	! - calculate diss production same in vector and cache optimization
37	! - move boundary condition for e and diss to boundary_conds
38	!
39	! 3545 2018-11-21 11:19:41Z gronemeier
40	! - Set rans_mode according to value of turbulence_closure
41	! - removed debug output
42	!
43	! 3430 2018-10-25 13:36:23Z maronga
44	! Added support for buildings in the dynamic SGS model
45	!
46	! 3398 2018-10-22 19:30:24Z knoop
47	! Refactored production_e and production_e_ij (removed code redundancy)
48	!
49	! 3386 2018-10-19 16:28:22Z gronemeier
50	! Renamed tcm_prognostic to tcm_prognostic_equations
51	!
52	! 3385 2018-10-19 14:52:29Z knoop
53	! Restructured loops and ifs in production_e to ease vectorization on GPUs
54	!
55	! 3300 2018-10-02 14:16:54Z gronemeier
56	! - removed global array wall_flags_0_global, hence reduced accuracy of l_wall
57	! calculation
58	! - removed maxloc call as this produced different results for different
59	! compiler options
60	!
61	! 3294 2018-10-01 02:37:10Z raasch
62	! changes concerning modularization of ocean option
63	!
64	! 3274 2018-09-24 15:42:55Z knoop
65	! Modularization of all bulk cloud physics code components
66	!
67	! 3245 2018-09-13 14:08:16Z knoop
68	! unused variables removed, shortest_distance has wp now
69	!
70	! 3183 2018-07-27 14:25:55Z suehring
71	! Rename variables and remove unused variable from USE statement
72	!
73	! 3182 2018-07-27 13:36:03Z suehring
74	! Use MOST for km only in RANS mode
75	!
76	! 3130 2018-07-16 11:08:55Z gronemeier
77	! - move boundary condition of km and kh to tcm_diffusivities
78	! - calculate km at boundaries according to MOST
79	! - move phi_m to surface_layer_fluxes_mod
80	!
81	! 3129 2018-07-16 07:45:13Z gronemeier
82	! - move limitation of diss to boundary_conds
83	! - move boundary conditions for e and diss to boundary_conds
84	! - consider non-default surfaces in tcm_diffusivities
85	! - use z_mo within surface layer instead of calculating it
86	! - resort output after case select -> reduced code duplication
87	! - when using rans_tke_e and 1d-model, do not use e1d, km1d and diss1d
88	!
89	! 3121 2018-07-11 18:46:49Z gronemeier
90	! - created the function phi_m
91	! - implemented km = u* * kappa * zp / phi_m in production_e_init for all
92	! surfaces
93	!
94	! 3120 2018-07-11 18:30:57Z gronemeier
95	! - changed tcm_diffusivities to tcm_diffusivities_default
96	! - created subroutine tcm_diffusivities that calls tcm_diffusivities_default
97	! and tcm_diffusivities_dynamic
98	!
99	! 3086 2018-06-25 09:08:04Z gronemeier
100	! bugfix: set rans_const_sigma(1) = 1.3
101	!
102	! 3083 2018-06-19 14:03:12Z gronemeier
103	! - set limits of diss at the end of prognostic equations
104	! - call production_e to calculate production term of diss
105	! - limit change of diss to -90% to +100%
106	! - remove factor 0.5 from diffusion_diss_ij
107	! - rename c_m into c_0, and c_h into c_4
108	! - add rans_const_c and rans_const_sigma as namelist parameters
109	! - add calculation of mixing length for profile output in case of rans_tke_e
110	! - changed format of annotations to comply with doxygen standards
111	! - calculate and save dissipation rate during rans_tke_l mode
112	! - set bc at vertical walls for e, diss, km, kh
113	! - bugfix: set l_wall = 0.0 within buildings
114	! - set l_wall at bottom and top boundary (rans-mode)
115	! - bugfix in production term for dissipation rate
116	! - bugfix in diffusion of dissipation rate
117	! - disable check for 1D model if rans_tke_e is used
118	! - bugfixes for initialization (rans-mode):
119	! - correction of dissipation-rate formula
120	! - calculate km based on l_wall
121	! - initialize diss if 1D model is not used
122	!
123	! 3045 2018-05-28 07:55:41Z Giersch
124	! Error message revised
125	!
126	! 3014 2018-05-09 08:42:38Z maronga
127	! Bugfix: nzb_do and nzt_do were not used for 3d data output
128	!
129	! 3004 2018-04-27 12:33:25Z Giersch
130	! Further allocation checks implemented
131	!
132	! 2938 2018-03-27 15:52:42Z suehring
133	! Further todo's
134	!
135	! 2936 2018-03-27 14:49:27Z gronemeier
136	! - defined l_grid only within this module
137	! - Moved l_wall definition from modules.f90
138	! - Get level of highest topography, used to limit upward distance calculation
139	! - Consider cyclic boundary conditions for mixing length calculation
140	! - Moved copy of wall_flags into subarray to subroutine
141	! - Implemented l_wall calculation in case of RANS simulation
142	! - Moved init of l_black to tcm_init_mixing_length
143	! - Moved init_mixing_length from init_grid.f90 and
144	! renamed it to tcm_init_mixing_length
145	!
146	! 2764 2018-01-22 09:25:36Z gronemeier
147	! Bugfix: remove duplicate SAVE statements
148	!
149	! 2746 2018-01-15 12:06:04Z suehring
150	! Move flag plant canopy to modules
151	!
152	! 2718 2018-01-02 08:49:38Z maronga
153	! Corrected "Former revisions" section
154	!
155	! 2701 2017-12-15 15:40:50Z suehring
156	! Changes from last commit documented
157	!
158	! 2698 2017-12-14 18:46:24Z suehring
159	! Bugfix in get_topography_top_index
160	!
161	! 2696 2017-12-14 17:12:51Z kanani
162	! Initial revision
163	!
164	!
165	!
166	!
167	! Authors:
168	! --------
169	! @author Tobias Gronemeier
170	! @author Hauke Wurps
171	!
172	! Description:
173	! ------------
174	!> This module contains the available turbulence closures for PALM.
175	!>
176	!>
177	!> @todo test initialization for all possibilities
178	!> add OpenMP directives whereever possible
179	!> @todo Check for random disturbances
180	!> @note <Enter notes on the module>
181	!------------------------------------------------------------------------------!
182	MODULE turbulence_closure_mod
183
184
185	USE arrays_3d, &
186	ONLY: diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3, &
187	e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m, &
188	te_m, tend, u, v, vpt, w
189
190	USE basic_constants_and_equations_mod, &
191	ONLY: g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
192
193	USE control_parameters, &
194	ONLY: constant_diffusion, dt_3d, e_init, humidity, &
195	initializing_actions, intermediate_timestep_count, &
196	intermediate_timestep_count_max, km_constant, &
197	les_dynamic, les_mw, ocean_mode, plant_canopy, prandtl_number, &
198	pt_reference, rans_mode, rans_tke_e, rans_tke_l, &
199	timestep_scheme, turbulence_closure, &
200	turbulent_inflow, use_upstream_for_tke, vpt_reference, &
201	ws_scheme_sca, current_timestep_number
202
203	USE advec_ws, &
204	ONLY: advec_s_ws
205
206	USE advec_s_bc_mod, &
207	ONLY: advec_s_bc
208
209	USE advec_s_pw_mod, &
210	ONLY: advec_s_pw
211
212	USE advec_s_up_mod, &
213	ONLY: advec_s_up
214
215	USE cpulog, &
216	ONLY: cpu_log, log_point, log_point_s
217
218	USE indices, &
219	ONLY: nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt, &
220	wall_flags_0
221
222	USE kinds
223
224	USE ocean_mod, &
225	ONLY: prho_reference
226
227	USE pegrid
228
229	USE plant_canopy_model_mod, &
230	ONLY: pcm_tendency
231
232	USE statistics, &
233	ONLY: hom, hom_sum, statistic_regions
234
235	USE user_actions_mod, &
236	ONLY: user_actions
237
238
239	IMPLICIT NONE
240
241
242	REAL(wp) :: c_0 !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
243	REAL(wp) :: c_1 !< model constant for RANS mode
244	REAL(wp) :: c_2 !< model constant for RANS mode
245	REAL(wp) :: c_3 !< model constant for RANS mode
246	REAL(wp) :: c_4 !< model constant for RANS mode
247	REAL(wp) :: l_max !< maximum length scale for Blackadar mixing length
248	REAL(wp) :: dsig_e = 1.0_wp !< factor to calculate Ke from Km (1/sigma_e)
249	REAL(wp) :: dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
250
251	REAL(wp), DIMENSION(0:4) :: rans_const_c = & !< model constants for RANS mode (namelist param)
252	(/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
253
254	REAL(wp), DIMENSION(2) :: rans_const_sigma = & !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
255	(/ 1.0_wp, 1.30_wp /)
256
257	REAL(wp), DIMENSION(:), ALLOCATABLE :: l_black !< mixing length according to Blackadar
258
259	REAL(wp), DIMENSION(:), ALLOCATABLE :: l_grid !< geometric mean of grid sizes dx, dy, dz
260
261	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: l_wall !< near-wall mixing length
262
263
264	PUBLIC c_0, rans_const_c, rans_const_sigma
265
266	!
267	!-- PALM interfaces:
268	!-- Input parameter checks to be done in check_parameters
269	INTERFACE tcm_check_parameters
270	MODULE PROCEDURE tcm_check_parameters
271	END INTERFACE tcm_check_parameters
272
273	!
274	!-- Data output checks for 2D/3D data to be done in check_parameters
275	INTERFACE tcm_check_data_output
276	MODULE PROCEDURE tcm_check_data_output
277	END INTERFACE tcm_check_data_output
278
279	!
280	!-- Definition of data output quantities
281	INTERFACE tcm_define_netcdf_grid
282	MODULE PROCEDURE tcm_define_netcdf_grid
283	END INTERFACE tcm_define_netcdf_grid
284
285	!
286	!-- Averaging of 3D data for output
287	INTERFACE tcm_3d_data_averaging
288	MODULE PROCEDURE tcm_3d_data_averaging
289	END INTERFACE tcm_3d_data_averaging
290
291	!
292	!-- Data output of 2D quantities
293	INTERFACE tcm_data_output_2d
294	MODULE PROCEDURE tcm_data_output_2d
295	END INTERFACE tcm_data_output_2d
296
297	!
298	!-- Data output of 3D data
299	INTERFACE tcm_data_output_3d
300	MODULE PROCEDURE tcm_data_output_3d
301	END INTERFACE tcm_data_output_3d
302
303	!
304	!-- Initialization actions
305	INTERFACE tcm_init
306	MODULE PROCEDURE tcm_init
307	END INTERFACE tcm_init
308
309	!
310	!-- Initialization of arrays
311	INTERFACE tcm_init_arrays
312	MODULE PROCEDURE tcm_init_arrays
313	END INTERFACE tcm_init_arrays
314
315	!
316	!-- Initialization of TKE production term
317	INTERFACE production_e_init
318	MODULE PROCEDURE production_e_init
319	END INTERFACE production_e_init
320
321	!
322	!-- Prognostic equations for TKE and TKE dissipation rate
323	INTERFACE tcm_prognostic_equations
324	MODULE PROCEDURE tcm_prognostic_equations
325	MODULE PROCEDURE tcm_prognostic_equations_ij
326	END INTERFACE tcm_prognostic_equations
327
328	!
329	!-- Production term for TKE
330	INTERFACE production_e
331	MODULE PROCEDURE production_e
332	MODULE PROCEDURE production_e_ij
333	END INTERFACE production_e
334
335	!
336	!-- Diffusion term for TKE
337	INTERFACE diffusion_e
338	MODULE PROCEDURE diffusion_e
339	MODULE PROCEDURE diffusion_e_ij
340	END INTERFACE diffusion_e
341
342	!
343	!-- Diffusion term for TKE dissipation rate
344	INTERFACE diffusion_diss
345	MODULE PROCEDURE diffusion_diss
346	MODULE PROCEDURE diffusion_diss_ij
347	END INTERFACE diffusion_diss
348
349	!
350	!-- Mixing length for LES case
351	INTERFACE mixing_length_les
352	MODULE PROCEDURE mixing_length_les
353	END INTERFACE mixing_length_les
354
355	!
356	!-- Mixing length for RANS case
357	INTERFACE mixing_length_rans
358	MODULE PROCEDURE mixing_length_rans
359	END INTERFACE mixing_length_rans
360
361	!
362	!-- Call tcm_diffusivities_default and tcm_diffusivities_dynamic
363	INTERFACE tcm_diffusivities
364	MODULE PROCEDURE tcm_diffusivities
365	END INTERFACE tcm_diffusivities
366
367	!
368	!-- Calculate diffusivities
369	INTERFACE tcm_diffusivities_default
370	MODULE PROCEDURE tcm_diffusivities_default
371	END INTERFACE tcm_diffusivities_default
372
373	!
374	!-- Calculate diffusivities according to dynamic sgs model
375	INTERFACE tcm_diffusivities_dynamic
376	MODULE PROCEDURE tcm_diffusivities_dynamic
377	END INTERFACE tcm_diffusivities_dynamic
378
379	!
380	!-- Box-filter method for dynamic sgs model
381	INTERFACE tcm_box_filter_2d
382	MODULE PROCEDURE tcm_box_filter_2d_single
383	MODULE PROCEDURE tcm_box_filter_2d_array
384	END INTERFACE tcm_box_filter_2d
385
386	!
387	!-- Swapping of time levels (required for prognostic variables)
388	INTERFACE tcm_swap_timelevel
389	MODULE PROCEDURE tcm_swap_timelevel
390	END INTERFACE tcm_swap_timelevel
391
392	SAVE
393
394	PRIVATE
395	!
396	!-- Add INTERFACES that must be available to other modules (alphabetical order)
397	PUBLIC production_e_init, tcm_3d_data_averaging, tcm_check_data_output, &
398	tcm_check_parameters, tcm_data_output_2d, tcm_data_output_3d, &
399	tcm_define_netcdf_grid, tcm_diffusivities, tcm_init, &
400	tcm_init_arrays, tcm_prognostic_equations, tcm_swap_timelevel
401
402
403	CONTAINS
404
405	!------------------------------------------------------------------------------!
406	! Description:
407	! ------------
408	!> Check parameters routine for turbulence closure module.
409	!------------------------------------------------------------------------------!
410	SUBROUTINE tcm_check_parameters
411
412	USE control_parameters, &
413	ONLY: message_string, turbulent_inflow, turbulent_outflow
414
415	IMPLICIT NONE
416
417	!
418	!-- Define which turbulence closure is going to be used
419	SELECT CASE ( TRIM( turbulence_closure ) )
420
421	CASE ( 'dynamic' )
422	les_dynamic = .TRUE.
423
424	CASE ( 'Moeng_Wyngaard' )
425	les_mw = .TRUE.
426
427	CASE ( 'TKE-l' )
428	rans_tke_l = .TRUE.
429	rans_mode = .TRUE.
430
431	CASE ( 'TKE-e' )
432	rans_tke_e = .TRUE.
433	rans_mode = .TRUE.
434
435	CASE DEFAULT
436	message_string = 'Unknown turbulence closure: ' // &
437	TRIM( turbulence_closure )
438	CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
439
440	END SELECT
441	!
442	!-- Set variables for RANS mode or LES mode
443	IF ( rans_mode ) THEN
444	!
445	!-- Assign values to constants for RANS mode
446	dsig_e = 1.0_wp / rans_const_sigma(1)
447	dsig_diss = 1.0_wp / rans_const_sigma(2)
448
449	c_0 = rans_const_c(0)
450	c_1 = rans_const_c(1)
451	c_2 = rans_const_c(2)
452	c_3 = rans_const_c(3) !> @todo clarify how to switch between different models
453	c_4 = rans_const_c(4)
454
455	IF ( turbulent_inflow .OR. turbulent_outflow ) THEN
456	message_string = 'turbulent inflow/outflow is not yet '// &
457	'implemented for RANS mode'
458	CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
459	ENDIF
460
461	message_string = 'RANS mode is still in development! ' // &
462	'&Not all features of PALM are yet compatible '// &
463	'with RANS mode. &Use at own risk!'
464	CALL message( 'tcm_check_parameters', 'PA0502', 0, 1, 0, 6, 0 )
465
466	ELSE
467	!
468	!-- LES mode
469	c_0 = 0.1_wp !according to Lilly (1967) and Deardorff (1980)
470
471	dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
472	!of K_e which is used in RANS mode
473
474	ENDIF
475
476	END SUBROUTINE tcm_check_parameters
477
478	!------------------------------------------------------------------------------!
479	! Description:
480	! ------------
481	!> Check data output.
482	!------------------------------------------------------------------------------!
483	SUBROUTINE tcm_check_data_output( var, unit )
484
485	IMPLICIT NONE
486
487	CHARACTER (LEN=*) :: unit !< unit of output variable
488	CHARACTER (LEN=*) :: var !< name of output variable
489
490
491	SELECT CASE ( TRIM( var ) )
492
493	CASE ( 'diss' )
494	unit = 'm2/s3'
495
496	CASE ( 'kh', 'km' )
497	unit = 'm2/s'
498
499	CASE DEFAULT
500	unit = 'illegal'
501
502	END SELECT
503
504	END SUBROUTINE tcm_check_data_output
505
506
507	!------------------------------------------------------------------------------!
508	! Description:
509	! ------------
510	!> Define appropriate grid for netcdf variables.
511	!> It is called out from subroutine netcdf.
512	!------------------------------------------------------------------------------!
513	SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
514
515	IMPLICIT NONE
516
517	CHARACTER (LEN=*), INTENT(OUT) :: grid_x !< x grid of output variable
518	CHARACTER (LEN=*), INTENT(OUT) :: grid_y !< y grid of output variable
519	CHARACTER (LEN=*), INTENT(OUT) :: grid_z !< z grid of output variable
520	CHARACTER (LEN=*), INTENT(IN) :: var !< name of output variable
521
522	LOGICAL, INTENT(OUT) :: found !< flag if output variable is found
523
524	found = .TRUE.
525
526	!
527	!-- Check for the grid
528	SELECT CASE ( TRIM( var ) )
529
530	CASE ( 'diss', 'diss_xy', 'diss_xz', 'diss_yz' )
531	grid_x = 'x'
532	grid_y = 'y'
533	grid_z = 'zu'
534
535	CASE ( 'kh', 'kh_xy', 'kh_xz', 'kh_yz' )
536	grid_x = 'x'
537	grid_y = 'y'
538	grid_z = 'zu'
539
540	CASE ( 'km', 'km_xy', 'km_xz', 'km_yz' )
541	grid_x = 'x'
542	grid_y = 'y'
543	grid_z = 'zu'
544
545	CASE DEFAULT
546	found = .FALSE.
547	grid_x = 'none'
548	grid_y = 'none'
549	grid_z = 'none'
550
551	END SELECT
552
553	END SUBROUTINE tcm_define_netcdf_grid
554
555
556	!------------------------------------------------------------------------------!
557	! Description:
558	! ------------
559	!> Average 3D data.
560	!------------------------------------------------------------------------------!
561	SUBROUTINE tcm_3d_data_averaging( mode, variable )
562
563
564	USE averaging, &
565	ONLY: diss_av, kh_av, km_av
566
567	USE control_parameters, &
568	ONLY: average_count_3d
569
570	IMPLICIT NONE
571
572	CHARACTER (LEN=*) :: mode !< flag defining mode 'allocate', 'sum' or 'average'
573	CHARACTER (LEN=*) :: variable !< name of variable
574
575	INTEGER(iwp) :: i !< loop index
576	INTEGER(iwp) :: j !< loop index
577	INTEGER(iwp) :: k !< loop index
578
579	IF ( mode == 'allocate' ) THEN
580
581	SELECT CASE ( TRIM( variable ) )
582
583	CASE ( 'diss' )
584	IF ( .NOT. ALLOCATED( diss_av ) ) THEN
585	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
586	ENDIF
587	diss_av = 0.0_wp
588
589	CASE ( 'kh' )
590	IF ( .NOT. ALLOCATED( kh_av ) ) THEN
591	ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
592	ENDIF
593	kh_av = 0.0_wp
594
595	CASE ( 'km' )
596	IF ( .NOT. ALLOCATED( km_av ) ) THEN
597	ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
598	ENDIF
599	km_av = 0.0_wp
600
601	CASE DEFAULT
602	CONTINUE
603
604	END SELECT
605
606	ELSEIF ( mode == 'sum' ) THEN
607
608	SELECT CASE ( TRIM( variable ) )
609
610	CASE ( 'diss' )
611	IF ( ALLOCATED( diss_av ) ) THEN
612	DO i = nxlg, nxrg
613	DO j = nysg, nyng
614	DO k = nzb, nzt+1
615	diss_av(k,j,i) = diss_av(k,j,i) + diss(k,j,i)
616	ENDDO
617	ENDDO
618	ENDDO
619	ENDIF
620
621	CASE ( 'kh' )
622	IF ( ALLOCATED( kh_av ) ) THEN
623	DO i = nxlg, nxrg
624	DO j = nysg, nyng
625	DO k = nzb, nzt+1
626	kh_av(k,j,i) = kh_av(k,j,i) + kh(k,j,i)
627	ENDDO
628	ENDDO
629	ENDDO
630	ENDIF
631
632	CASE ( 'km' )
633	IF ( ALLOCATED( km_av ) ) THEN
634	DO i = nxlg, nxrg
635	DO j = nysg, nyng
636	DO k = nzb, nzt+1
637	km_av(k,j,i) = km_av(k,j,i) + km(k,j,i)
638	ENDDO
639	ENDDO
640	ENDDO
641	ENDIF
642
643	CASE DEFAULT
644	CONTINUE
645
646	END SELECT
647
648	ELSEIF ( mode == 'average' ) THEN
649
650	SELECT CASE ( TRIM( variable ) )
651
652	CASE ( 'diss' )
653	IF ( ALLOCATED( diss_av ) ) THEN
654	DO i = nxlg, nxrg
655	DO j = nysg, nyng
656	DO k = nzb, nzt+1
657	diss_av(k,j,i) = diss_av(k,j,i) &
658	/ REAL( average_count_3d, KIND=wp )
659	ENDDO
660	ENDDO
661	ENDDO
662	ENDIF
663
664	CASE ( 'kh' )
665	IF ( ALLOCATED( kh_av ) ) THEN
666	DO i = nxlg, nxrg
667	DO j = nysg, nyng
668	DO k = nzb, nzt+1
669	kh_av(k,j,i) = kh_av(k,j,i) &
670	/ REAL( average_count_3d, KIND=wp )
671	ENDDO
672	ENDDO
673	ENDDO
674	ENDIF
675
676	CASE ( 'km' )
677	IF ( ALLOCATED( km_av ) ) THEN
678	DO i = nxlg, nxrg
679	DO j = nysg, nyng
680	DO k = nzb, nzt+1
681	km_av(k,j,i) = km_av(k,j,i) &
682	/ REAL( average_count_3d, KIND=wp )
683	ENDDO
684	ENDDO
685	ENDDO
686	ENDIF
687
688	END SELECT
689
690	ENDIF
691
692	END SUBROUTINE tcm_3d_data_averaging
693
694
695	!------------------------------------------------------------------------------!
696	! Description:
697	! ------------
698	!> Define 2D output variables.
699	!------------------------------------------------------------------------------!
700	SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, &
701	nzb_do, nzt_do )
702
703	USE averaging, &
704	ONLY: diss_av, kh_av, km_av
705
706	IMPLICIT NONE
707
708	CHARACTER (LEN=*) :: grid !< name of vertical grid
709	CHARACTER (LEN=*) :: mode !< either 'xy', 'xz' or 'yz'
710	CHARACTER (LEN=*) :: variable !< name of variable
711
712	INTEGER(iwp) :: av !< flag for (non-)average output
713	INTEGER(iwp) :: flag_nr !< number of masking flag
714	INTEGER(iwp) :: i !< loop index
715	INTEGER(iwp) :: j !< loop index
716	INTEGER(iwp) :: k !< loop index
717	INTEGER(iwp) :: nzb_do !< vertical output index (bottom)
718	INTEGER(iwp) :: nzt_do !< vertical output index (top)
719
720	LOGICAL :: found !< flag if output variable is found
721	LOGICAL :: resorted !< flag if output is already resorted
722
723	REAL(wp) :: fill_value = -9999.0_wp !< value for the _FillValue attribute
724
725	REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
726	!< array to which output data is resorted to
727
728	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
729
730	found = .TRUE.
731	resorted = .FALSE.
732	!
733	!-- Set masking flag for topography for not resorted arrays
734	flag_nr = 0
735
736	SELECT CASE ( TRIM( variable ) )
737
738	CASE ( 'diss_xy', 'diss_xz', 'diss_yz' )
739	IF ( av == 0 ) THEN
740	to_be_resorted => diss
741	ELSE
742	IF ( .NOT. ALLOCATED( diss_av ) ) THEN
743	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
744	diss_av = REAL( fill_value, KIND = wp )
745	ENDIF
746	to_be_resorted => diss_av
747	ENDIF
748	IF ( mode == 'xy' ) grid = 'zu'
749
750	CASE ( 'kh_xy', 'kh_xz', 'kh_yz' )
751	IF ( av == 0 ) THEN
752	to_be_resorted => kh
753	ELSE
754	IF ( .NOT. ALLOCATED( kh_av ) ) THEN
755	ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
756	kh_av = REAL( fill_value, KIND = wp )
757	ENDIF
758	to_be_resorted => kh_av
759	ENDIF
760	IF ( mode == 'xy' ) grid = 'zu'
761
762	CASE ( 'km_xy', 'km_xz', 'km_yz' )
763	IF ( av == 0 ) THEN
764	to_be_resorted => km
765	ELSE
766	IF ( .NOT. ALLOCATED( km_av ) ) THEN
767	ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
768	km_av = REAL( fill_value, KIND = wp )
769	ENDIF
770	to_be_resorted => km_av
771	ENDIF
772	IF ( mode == 'xy' ) grid = 'zu'
773
774	CASE DEFAULT
775	found = .FALSE.
776	grid = 'none'
777
778	END SELECT
779
780	IF ( found .AND. .NOT. resorted ) THEN
781	DO i = nxl, nxr
782	DO j = nys, nyn
783	DO k = nzb_do, nzt_do
784	local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), &
785	REAL( fill_value, KIND = wp ), &
786	BTEST( wall_flags_0(k,j,i), flag_nr ) )
787	ENDDO
788	ENDDO
789	ENDDO
790	ENDIF
791
792	END SUBROUTINE tcm_data_output_2d
793
794
795	!------------------------------------------------------------------------------!
796	! Description:
797	! ------------
798	!> Define 3D output variables.
799	!------------------------------------------------------------------------------!
800	SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
801
802
803	USE averaging, &
804	ONLY: diss_av, kh_av, km_av
805
806	IMPLICIT NONE
807
808	CHARACTER (LEN=*) :: variable !< name of variable
809
810	INTEGER(iwp) :: av !< flag for (non-)average output
811	INTEGER(iwp) :: flag_nr !< number of masking flag
812	INTEGER(iwp) :: i !< loop index
813	INTEGER(iwp) :: j !< loop index
814	INTEGER(iwp) :: k !< loop index
815	INTEGER(iwp) :: nzb_do !< lower limit of the data output (usually 0)
816	INTEGER(iwp) :: nzt_do !< vertical upper limit of the data output (usually nz_do3d)
817
818	LOGICAL :: found !< flag if output variable is found
819	LOGICAL :: resorted !< flag if output is already resorted
820
821	REAL(wp) :: fill_value = -9999.0_wp !< value for the _FillValue attribute
822
823	REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) :: local_pf !< local
824	!< array to which output data is resorted to
825
826	REAL(wp), DIMENSION(:,:,:), POINTER :: to_be_resorted !< points to selected output variable
827
828	found = .TRUE.
829	resorted = .FALSE.
830	!
831	!-- Set masking flag for topography for not resorted arrays
832	flag_nr = 0
833
834	SELECT CASE ( TRIM( variable ) )
835
836	CASE ( 'diss' )
837	IF ( av == 0 ) THEN
838	to_be_resorted => diss
839	ELSE
840	IF ( .NOT. ALLOCATED( diss_av ) ) THEN
841	ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
842	diss_av = REAL( fill_value, KIND = wp )
843	ENDIF
844	to_be_resorted => diss_av
845	ENDIF
846
847	CASE ( 'kh' )
848	IF ( av == 0 ) THEN
849	to_be_resorted => kh
850	ELSE
851	IF ( .NOT. ALLOCATED( kh_av ) ) THEN
852	ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
853	kh_av = REAL( fill_value, KIND = wp )
854	ENDIF
855	to_be_resorted => kh_av
856	ENDIF
857
858	CASE ( 'km' )
859	IF ( av == 0 ) THEN
860	to_be_resorted => km
861	ELSE
862	IF ( .NOT. ALLOCATED( km_av ) ) THEN
863	ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
864	km_av = REAL( fill_value, KIND = wp )
865	ENDIF
866	to_be_resorted => km_av
867	ENDIF
868
869	CASE DEFAULT
870	found = .FALSE.
871
872	END SELECT
873
874
875	IF ( found .AND. .NOT. resorted ) THEN
876	DO i = nxl, nxr
877	DO j = nys, nyn
878	DO k = nzb_do, nzt_do
879	local_pf(i,j,k) = MERGE( &
880	to_be_resorted(k,j,i), &
881	REAL( fill_value, KIND = wp ), &
882	BTEST( wall_flags_0(k,j,i), flag_nr ) )
883	ENDDO
884	ENDDO
885	ENDDO
886	resorted = .TRUE.
887	ENDIF
888
889	END SUBROUTINE tcm_data_output_3d
890
891
892	!------------------------------------------------------------------------------!
893	! Description:
894	! ------------
895	!> Allocate arrays and assign pointers.
896	!------------------------------------------------------------------------------!
897	SUBROUTINE tcm_init_arrays
898
899	USE bulk_cloud_model_mod, &
900	ONLY: collision_turbulence
901
902	USE particle_attributes, &
903	ONLY: use_sgs_for_particles, wang_kernel
904
905	USE pmc_interface, &
906	ONLY: nested_run
907
908	IMPLICIT NONE
909
910	ALLOCATE( kh(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
911	ALLOCATE( km(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
912
913	ALLOCATE( e_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
914	ALLOCATE( e_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
915	ALLOCATE( e_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
916
917	!
918	!-- Allocate arrays required for dissipation.
919	!-- Please note, if it is a nested run, arrays need to be allocated even if
920	!-- they do not necessarily need to be transferred, which is attributed to
921	!-- the design of the model coupler which allocates memory for each variable.
922	IF ( rans_mode .OR. use_sgs_for_particles .OR. wang_kernel .OR. &
923	collision_turbulence .OR. nested_run ) THEN
924
925	ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
926	IF ( rans_tke_e .OR. nested_run ) THEN
927	ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
928	ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
929	ENDIF
930
931	ENDIF
932
933	!
934	!-- Initial assignment of pointers
935	e => e_1; e_p => e_2; te_m => e_3
936
937	IF ( rans_mode .OR. use_sgs_for_particles .OR. &
938	wang_kernel .OR. collision_turbulence .OR. nested_run ) THEN
939	diss => diss_1
940	IF ( rans_tke_e .OR. nested_run ) THEN
941	diss_p => diss_2; tdiss_m => diss_3
942	ENDIF
943	ENDIF
944
945	END SUBROUTINE tcm_init_arrays
946
947
948	!------------------------------------------------------------------------------!
949	! Description:
950	! ------------
951	!> Initialization of turbulence closure module.
952	!------------------------------------------------------------------------------!
953	SUBROUTINE tcm_init
954
955	USE control_parameters, &
956	ONLY: bc_dirichlet_l, complex_terrain, topography
957
958	USE model_1d_mod, &
959	ONLY: e1d, kh1d, km1d
960
961	USE surface_mod, &
962	ONLY: get_topography_top_index_ji
963
964	IMPLICIT NONE
965
966	INTEGER(iwp) :: i !< loop index
967	INTEGER(iwp) :: j !< loop index
968	INTEGER(iwp) :: k !< loop index
969	INTEGER(iwp) :: nz_s_shift !< lower shift index for scalars
970	INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
971
972	!
973	!-- Initialize mixing length
974	CALL tcm_init_mixing_length
975
976	!
977	!-- Actions for initial runs
978	IF ( TRIM( initializing_actions ) /= 'read_restart_data' .AND. &
979	TRIM( initializing_actions ) /= 'cyclic_fill' ) THEN
980
981	IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 ) THEN
982
983	IF ( .NOT. rans_tke_e ) THEN
984	!
985	!-- Transfer initial profiles to the arrays of the 3D model
986	DO i = nxlg, nxrg
987	DO j = nysg, nyng
988	e(:,j,i) = e1d
989	kh(:,j,i) = kh1d
990	km(:,j,i) = km1d
991	ENDDO
992	ENDDO
993
994	IF ( constant_diffusion ) THEN
995	e = 0.0_wp
996	ENDIF
997
998	ELSE
999	!
1000	!-- In case of TKE-e closure in RANS mode, do not use e, diss, and km
1001	!-- profiles from 1D model. Instead, initialize with constant profiles
1002	IF ( constant_diffusion ) THEN
1003	km = km_constant
1004	kh = km / prandtl_number
1005	e = 0.0_wp
1006	ELSEIF ( e_init > 0.0_wp ) THEN
1007	DO i = nxlg, nxrg
1008	DO j = nysg, nyng
1009	DO k = nzb+1, nzt
1010	km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1011	ENDDO
1012	ENDDO
1013	ENDDO
1014	km(nzb,:,:) = km(nzb+1,:,:)
1015	km(nzt+1,:,:) = km(nzt,:,:)
1016	kh = km / prandtl_number
1017	e = e_init
1018	ELSE
1019	IF ( .NOT. ocean_mode ) THEN
1020	kh = 0.01_wp ! there must exist an initial diffusion, because
1021	km = 0.01_wp ! otherwise no TKE would be produced by the
1022	! production terms, as long as not yet
1023	! e = (u/cm)*2 at k=nzb+1
1024	ELSE
1025	kh = 0.00001_wp
1026	km = 0.00001_wp
1027	ENDIF
1028	e = 0.0_wp
1029	ENDIF
1030
1031	DO i = nxlg, nxrg
1032	DO j = nysg, nyng
1033	DO k = nzb+1, nzt
1034	diss(k,j,i) = c_0*4 e(k,j,i)**2 / km(k,j,i)
1035	ENDDO
1036	ENDDO
1037	ENDDO
1038	diss(nzb,:,:) = diss(nzb+1,:,:)
1039	diss(nzt+1,:,:) = diss(nzt,:,:)
1040
1041	ENDIF
1042
1043	ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 .OR. &
1044	INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN
1045
1046	IF ( constant_diffusion ) THEN
1047	km = km_constant
1048	kh = km / prandtl_number
1049	e = 0.0_wp
1050	ELSEIF ( e_init > 0.0_wp ) THEN
1051	DO i = nxlg, nxrg
1052	DO j = nysg, nyng
1053	DO k = nzb+1, nzt
1054	km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
1055	ENDDO
1056	ENDDO
1057	ENDDO
1058	km(nzb,:,:) = km(nzb+1,:,:)
1059	km(nzt+1,:,:) = km(nzt,:,:)
1060	kh = km / prandtl_number
1061	e = e_init
1062	ELSE
1063	IF ( .NOT. ocean_mode ) THEN
1064	kh = 0.01_wp ! there must exist an initial diffusion, because
1065	km = 0.01_wp ! otherwise no TKE would be produced by the
1066	! production terms, as long as not yet
1067	! e = (u/cm)*2 at k=nzb+1
1068	ELSE
1069	kh = 0.00001_wp
1070	km = 0.00001_wp
1071	ENDIF
1072	e = 0.0_wp
1073	ENDIF
1074
1075	IF ( rans_tke_e ) THEN
1076	DO i = nxlg, nxrg
1077	DO j = nysg, nyng
1078	DO k = nzb+1, nzt
1079	diss(k,j,i) = c_0*4 e(k,j,i)**2 / km(k,j,i)
1080	ENDDO
1081	ENDDO
1082	ENDDO
1083	diss(nzb,:,:) = diss(nzb+1,:,:)
1084	diss(nzt+1,:,:) = diss(nzt,:,:)
1085	ENDIF
1086
1087	ENDIF
1088	!
1089	!-- Store initial profiles for output purposes etc.
1090	hom(:,1,23,:) = SPREAD( km(:,nys,nxl), 2, statistic_regions+1 )
1091	hom(:,1,24,:) = SPREAD( kh(:,nys,nxl), 2, statistic_regions+1 )
1092	!
1093	!-- Initialize old and new time levels.
1094	te_m = 0.0_wp
1095	e_p = e
1096	IF ( rans_tke_e ) THEN
1097	tdiss_m = 0.0_wp
1098	diss_p = diss
1099	ENDIF
1100
1101	ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data' .OR. &
1102	TRIM( initializing_actions ) == 'cyclic_fill' ) &
1103	THEN
1104
1105	!
1106	!-- In case of complex terrain and cyclic fill method as initialization,
1107	!-- shift initial data in the vertical direction for each point in the
1108	!-- x-y-plane depending on local surface height
1109	IF ( complex_terrain .AND. &
1110	TRIM( initializing_actions ) == 'cyclic_fill' ) THEN
1111	DO i = nxlg, nxrg
1112	DO j = nysg, nyng
1113	nz_s_shift = get_topography_top_index_ji( j, i, 's' )
1114
1115	e(nz_s_shift:nzt+1,j,i) = e(0:nzt+1-nz_s_shift,j,i)
1116	km(nz_s_shift:nzt+1,j,i) = km(0:nzt+1-nz_s_shift,j,i)
1117	kh(nz_s_shift:nzt+1,j,i) = kh(0:nzt+1-nz_s_shift,j,i)
1118	ENDDO
1119	ENDDO
1120	IF ( rans_tke_e ) THEN
1121	DO i = nxlg, nxrg
1122	DO j = nysg, nyng
1123	nz_s_shift = get_topography_top_index_ji( j, i, 's' )
1124
1125	diss(nz_s_shift:nzt+1,j,i) = diss(0:nzt+1-nz_s_shift,j,i)
1126	ENDDO
1127	ENDDO
1128	ENDIF
1129	ENDIF
1130
1131	!
1132	!-- Initialization of the turbulence recycling method
1133	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
1134	turbulent_inflow ) THEN
1135	mean_inflow_profiles(:,5) = hom_sum(:,8,0) ! e
1136	!
1137	!-- In case of complex terrain, determine vertical displacement at inflow
1138	!-- boundary and adjust mean inflow profiles
1139	IF ( complex_terrain ) THEN
1140	IF ( nxlg <= 0 .AND. nxrg >= 0 .AND. &
1141	nysg <= 0 .AND. nyng >= 0 ) THEN
1142	nz_s_shift_l = get_topography_top_index_ji( 0, 0, 's' )
1143	ELSE
1144	nz_s_shift_l = 0
1145	ENDIF
1146	#if defined( __parallel )
1147	CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, &
1148	MPI_MAX, comm2d, ierr)
1149	#else
1150	nz_s_shift = nz_s_shift_l
1151	#endif
1152	mean_inflow_profiles(nz_s_shift:nzt+1,5) = &
1153	hom_sum(0:nzt+1-nz_s_shift,8,0) ! e
1154	ENDIF
1155	!
1156	!-- Use these mean profiles at the inflow (provided that Dirichlet
1157	!-- conditions are used)
1158	IF ( bc_dirichlet_l ) THEN
1159	DO j = nysg, nyng
1160	DO k = nzb, nzt+1
1161	e(k,j,nxlg:-1) = mean_inflow_profiles(k,5)
1162	ENDDO
1163	ENDDO
1164	ENDIF
1165	ENDIF
1166	!
1167	!-- Inside buildings set TKE back to zero
1168	IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. &
1169	topography /= 'flat' ) THEN
1170	!
1171	!-- Inside buildings set TKE back to zero.
1172	!-- Other scalars (km, kh,...) are ignored at present,
1173	!-- maybe revise later.
1174	DO i = nxlg, nxrg
1175	DO j = nysg, nyng
1176	DO k = nzb, nzt
1177	e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, &
1178	BTEST( wall_flags_0(k,j,i), 0 ) )
1179	ENDDO
1180	ENDDO
1181	ENDDO
1182
1183	IF ( rans_tke_e ) THEN
1184	DO i = nxlg, nxrg
1185	DO j = nysg, nyng
1186	DO k = nzb, nzt
1187	diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp, &
1188	BTEST( wall_flags_0(k,j,i), 0 ) )
1189	ENDDO
1190	ENDDO
1191	ENDDO
1192	ENDIF
1193	ENDIF
1194	!
1195	!-- Initialize new time levels (only done in order to set boundary values
1196	!-- including ghost points)
1197	e_p = e
1198	!
1199	!-- Allthough tendency arrays are set in prognostic_equations, they have
1200	!-- to be predefined here because there they are used (but multiplied with 0)
1201	!-- before they are set.
1202	te_m = 0.0_wp
1203
1204	IF ( rans_tke_e ) THEN
1205	diss_p = diss
1206	tdiss_m = 0.0_wp
1207	ENDIF
1208
1209	ENDIF
1210
1211	END SUBROUTINE tcm_init
1212
1213
1214	! Description:
1215	! -----------------------------------------------------------------------------!
1216	!> Pre-computation of grid-dependent and near-wall mixing length.
1217	!> @todo consider walls in horizontal direction at a distance further than a
1218	!> single grid point (RANS mode)
1219	!------------------------------------------------------------------------------!
1220	SUBROUTINE tcm_init_mixing_length
1221
1222	USE arrays_3d, &
1223	ONLY: dzw, ug, vg, zu, zw
1224
1225	USE control_parameters, &
1226	ONLY: bc_lr_cyc, bc_ns_cyc, f, message_string, wall_adjustment_factor
1227
1228	USE grid_variables, &
1229	ONLY: dx, dy
1230
1231	USE indices, &
1232	ONLY: nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, &
1233	nzt, wall_flags_0
1234
1235	USE kinds
1236
1237
1238	IMPLICIT NONE
1239
1240	INTEGER(iwp) :: dist_dx !< found distance devided by dx
1241	INTEGER(iwp) :: i !< index variable along x
1242	INTEGER(iwp) :: ii !< index variable along x
1243	INTEGER(iwp) :: j !< index variable along y
1244	INTEGER(iwp) :: k !< index variable along z
1245	INTEGER(iwp) :: k_max_topo = 0 !< index of maximum topography height
1246	INTEGER(iwp) :: kk !< index variable along z
1247	INTEGER(iwp) :: rad_i !< search radius in grid points along x
1248	INTEGER(iwp) :: rad_i_l !< possible search radius to the left
1249	INTEGER(iwp) :: rad_i_r !< possible search radius to the right
1250	INTEGER(iwp) :: rad_j !< search radius in grid points along y
1251	INTEGER(iwp) :: rad_j_n !< possible search radius to north
1252	INTEGER(iwp) :: rad_j_s !< possible search radius to south
1253	INTEGER(iwp) :: rad_k !< search radius in grid points along z
1254	INTEGER(iwp) :: rad_k_b !< search radius in grid points along negative z
1255	INTEGER(iwp) :: rad_k_t !< search radius in grid points along positive z
1256
1257	INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
1258
1259	INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
1260
1261	INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE :: wall_flags_dummy !< dummy array required for MPI_ALLREDUCE command
1262
1263	REAL(wp) :: radius !< search radius in meter
1264
1265	ALLOCATE( l_grid(1:nzt) )
1266	ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1267	!
1268	!-- Initialize the mixing length in case of an LES-simulation
1269	IF ( .NOT. rans_mode ) THEN
1270	!
1271	!-- Compute the grid-dependent mixing length.
1272	DO k = 1, nzt
1273	l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp
1274	ENDDO
1275	!
1276	!-- Initialize near-wall mixing length l_wall only in the vertical direction
1277	!-- for the moment, multiplication with wall_adjustment_factor further below
1278	l_wall(nzb,:,:) = l_grid(1)
1279	DO k = nzb+1, nzt
1280	l_wall(k,:,:) = l_grid(k)
1281	ENDDO
1282	l_wall(nzt+1,:,:) = l_grid(nzt)
1283
1284	DO k = 1, nzt
1285	IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR. &
1286	l_grid(k) > 1.5_wp * dy * wall_adjustment_factor ) THEN
1287	WRITE( message_string, * ) 'grid anisotropy exceeds ', &
1288	'threshold given by only local', &
1289	' &horizontal reduction of near_wall ', &
1290	'mixing length l_wall', &
1291	' &starting from height level k = ', k, &
1292	'.'
1293	CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
1294	EXIT
1295	ENDIF
1296	ENDDO
1297	!
1298	!-- In case of topography: limit near-wall mixing length l_wall further:
1299	!-- Go through all points of the subdomain one by one and look for the closest
1300	!-- surface.
1301	!-- Is this correct in the ocean case?
1302	DO i = nxl, nxr
1303	DO j = nys, nyn
1304	DO k = nzb+1, nzt
1305	!
1306	!-- Check if current gridpoint belongs to the atmosphere
1307	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1308	!
1309	!-- Check for neighbouring grid-points.
1310	!-- Vertical distance, down
1311	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) ) &
1312	l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
1313	!
1314	!-- Vertical distance, up
1315	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) ) &
1316	l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
1317	!
1318	!-- y-distance
1319	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 ) .OR. &
1320	.NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) ) &
1321	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
1322	!
1323	!-- x-distance
1324	IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 ) .OR. &
1325	.NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) ) &
1326	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
1327	!
1328	!-- yz-distance (vertical edges, down)
1329	IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 ) .OR. &
1330	.NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 ) ) &
1331	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1332	SQRT( 0.25_wp * dy**2 + &
1333	( zu(k) - zw(k-1) )**2 ) )
1334	!
1335	!-- yz-distance (vertical edges, up)
1336	IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 ) .OR. &
1337	.NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 ) ) &
1338	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1339	SQRT( 0.25_wp * dy**2 + &
1340	( zw(k) - zu(k) )**2 ) )
1341	!
1342	!-- xz-distance (vertical edges, down)
1343	IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 ) .OR. &
1344	.NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 ) ) &
1345	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1346	SQRT( 0.25_wp * dx**2 + &
1347	( zu(k) - zw(k-1) )**2 ) )
1348	!
1349	!-- xz-distance (vertical edges, up)
1350	IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 ) .OR. &
1351	.NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 ) ) &
1352	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1353	SQRT( 0.25_wp * dx**2 + &
1354	( zw(k) - zu(k) )**2 ) )
1355	!
1356	!-- xy-distance (horizontal edges)
1357	IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 ) .OR. &
1358	.NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 ) .OR. &
1359	.NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 ) .OR. &
1360	.NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) ) &
1361	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1362	SQRT( 0.25_wp * ( dx2 + dy2 ) ) )
1363	!
1364	!-- xyz distance (vertical and horizontal edges, down)
1365	IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 ) .OR. &
1366	.NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 ) .OR. &
1367	.NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 ) .OR. &
1368	.NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) ) &
1369	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1370	SQRT( 0.25_wp * ( dx2 + dy2 ) &
1371	+ ( zu(k) - zw(k-1) )**2 ) )
1372	!
1373	!-- xyz distance (vertical and horizontal edges, up)
1374	IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 ) .OR. &
1375	.NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 ) .OR. &
1376	.NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 ) .OR. &
1377	.NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) ) &
1378	l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), &
1379	SQRT( 0.25_wp * ( dx2 + dy2 ) &
1380	+ ( zw(k) - zu(k) )**2 ) )
1381
1382	ENDIF
1383	ENDDO
1384	ENDDO
1385	ENDDO
1386
1387	ELSE
1388	!
1389	!-- Initialize the mixing length in case of a RANS simulation
1390	ALLOCATE( l_black(nzb:nzt+1) )
1391
1392	!
1393	!-- Calculate mixing length according to Blackadar (1962)
1394	IF ( f /= 0.0_wp ) THEN
1395	l_max = 2.7E-4_wp * SQRT( ug(nzt+1)2 + vg(nzt+1)2 ) / &
1396	ABS( f ) + 1.0E-10_wp
1397	ELSE
1398	l_max = 30.0_wp
1399	ENDIF
1400
1401	DO k = nzb, nzt
1402	l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max )
1403	ENDDO
1404
1405	l_black(nzt+1) = l_black(nzt)
1406
1407	!
1408	!-- Get height level of highest topography within local subdomain
1409	DO i = nxlg, nxrg
1410	DO j = nysg, nyng
1411	DO k = nzb+1, nzt-1
1412	IF ( .NOT. BTEST( wall_flags_0(k,j,i), 0 ) .AND. &
1413	k > k_max_topo ) &
1414	k_max_topo = k
1415	ENDDO
1416	ENDDO
1417	ENDDO
1418
1419	l_wall(nzb,:,:) = l_black(nzb)
1420	l_wall(nzt+1,:,:) = l_black(nzt+1)
1421	!
1422	!-- Limit mixing length to either nearest wall or Blackadar mixing length.
1423	!-- For that, analyze each grid point (i/j/k) ("analysed grid point") and
1424	!-- search within its vicinity for the shortest distance to a wall by cal-
1425	!-- culating the distance between the analysed grid point and the "viewed
1426	!-- grid point" if it contains a wall (belongs to topography).
1427	DO k = nzb+1, nzt
1428
1429	radius = l_black(k) ! radius within walls are searched
1430	!
1431	!-- Set l_wall to its default maximum value (l_back)
1432	l_wall(k,:,:) = radius
1433
1434	!
1435	!-- Compute search radius as number of grid points in all directions
1436	rad_i = CEILING( radius / dx )
1437	rad_j = CEILING( radius / dy )
1438
1439	DO kk = 0, nzt-k
1440	rad_k_t = kk
1441	!
1442	!-- Limit upward search radius to height of maximum topography
1443	IF ( zu(k+kk)-zu(k) >= radius .OR. k+kk >= k_max_topo ) EXIT
1444	ENDDO
1445
1446	DO kk = 0, k
1447	rad_k_b = kk
1448	IF ( zu(k)-zu(k-kk) >= radius ) EXIT
1449	ENDDO
1450
1451	!
1452	!-- Get maximum vertical radius; necessary for defining arrays
1453	rad_k = MAX( rad_k_b, rad_k_t )
1454	!
1455	!-- When analysed grid point lies above maximum topography, set search
1456	!-- radius to 0 if the distance between the analysed grid point and max
1457	!-- topography height is larger than the maximum search radius
1458	IF ( zu(k-rad_k_b) > zu(k_max_topo) ) rad_k_b = 0
1459	!
1460	!-- Search within vicinity only if the vertical search radius is >0
1461	IF ( rad_k_b /= 0 .OR. rad_k_t /= 0 ) THEN
1462
1463	!> @note shape of vicinity is larger in z direction
1464	!> Shape of vicinity is two grid points larger than actual search
1465	!> radius in vertical direction. The first and last grid point is
1466	!> always set to 1 to asure correct detection of topography. See
1467	!> function "shortest_distance" for details.
1468	!> 2018-03-16, gronemeier
1469	ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
1470	ALLOCATE( vic_yz(0:rad_k+1,0:rad_j) )
1471
1472	vicinity = 1
1473
1474	DO i = nxl, nxr
1475	DO j = nys, nyn
1476	!
1477	!-- Start search only if (i/j/k) belongs to atmosphere
1478	IF ( BTEST( wall_flags_0(k,j,i), 0 ) ) THEN
1479	!
1480	!-- Reset topography within vicinity
1481	vicinity(-rad_k:rad_k,:,:) = 0
1482	!
1483	!-- Copy area surrounding analysed grid point into vicinity.
1484	!-- First, limit size of data copied to vicinity by the domain
1485	!-- border
1486	!> @note limit copied area to 1 grid point in hor. dir.
1487	!> Ignore walls in horizontal direction which are
1488	!> further away than a single grid point. This allows to
1489	!> only search within local subdomain without the need
1490	!> of global topography information.
1491	!> The error made by this assumption are acceptable at
1492	!> the moment.
1493	!> 2018-10-01, gronemeier
1494	rad_i_l = MIN( 1, rad_i, i )
1495	rad_i_r = MIN( 1, rad_i, nx-i )
1496
1497	rad_j_s = MIN( 1, rad_j, j )
1498	rad_j_n = MIN( 1, rad_j, ny-j )
1499
1500	CALL copy_into_vicinity( k, j, i, &
1501	-rad_k_b, rad_k_t, &
1502	-rad_j_s, rad_j_n, &
1503	-rad_i_l, rad_i_r )
1504	!> @note in case of cyclic boundaries, those parts of the
1505	!> topography which lies beyond the domain borders but
1506	!> still within the search radius still needs to be
1507	!> copied into 'vicinity'. As the effective search
1508	!> radius is limited to 1 at the moment, no further
1509	!> copying is needed. Old implementation (prior to
1510	!> 2018-10-01) had this covered but used a global array.
1511	!> 2018-10-01, gronemeier
1512
1513	!
1514	!-- Search for walls only if there is any within vicinity
1515	IF ( MAXVAL( vicinity(-rad_k:rad_k,:,:) ) /= 0 ) THEN
1516	!
1517	!-- Search within first half (positive x)
1518	dist_dx = rad_i
1519	DO ii = 0, dist_dx
1520	!
1521	!-- Search along vertical direction only if below
1522	!-- maximum topography
1523	IF ( rad_k_t > 0 ) THEN
1524	!
1525	!-- Search for walls within octant (+++)
1526	vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1527	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1528	shortest_distance( vic_yz, .TRUE., ii ) )
1529	!
1530	!-- Search for walls within octant (+-+)
1531	!-- Switch order of array so that the analysed grid
1532	!-- point is always located at (0/0) (required by
1533	!-- shortest_distance").
1534	vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1535	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1536	shortest_distance( vic_yz, .TRUE., ii ) )
1537
1538	ENDIF
1539	!
1540	!-- Search for walls within octant (+--)
1541	vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1542	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1543	shortest_distance( vic_yz, .FALSE., ii ) )
1544	!
1545	!-- Search for walls within octant (++-)
1546	vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1547	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1548	shortest_distance( vic_yz, .FALSE., ii ) )
1549	!
1550	!-- Reduce search along x by already found distance
1551	dist_dx = CEILING( l_wall(k,j,i) / dx )
1552
1553	ENDDO
1554	!
1555	!- Search within second half (negative x)
1556	DO ii = 0, -dist_dx, -1
1557	!
1558	!-- Search along vertical direction only if below
1559	!-- maximum topography
1560	IF ( rad_k_t > 0 ) THEN
1561	!
1562	!-- Search for walls within octant (-++)
1563	vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
1564	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1565	shortest_distance( vic_yz, .TRUE., -ii ) )
1566	!
1567	!-- Search for walls within octant (--+)
1568	!-- Switch order of array so that the analysed grid
1569	!-- point is always located at (0/0) (required by
1570	!-- shortest_distance").
1571	vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
1572	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1573	shortest_distance( vic_yz, .TRUE., -ii ) )
1574
1575	ENDIF
1576	!
1577	!-- Search for walls within octant (---)
1578	vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
1579	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1580	shortest_distance( vic_yz, .FALSE., -ii ) )
1581	!
1582	!-- Search for walls within octant (-+-)
1583	vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
1584	l_wall(k,j,i) = MIN( l_wall(k,j,i), &
1585	shortest_distance( vic_yz, .FALSE., -ii ) )
1586	!
1587	!-- Reduce search along x by already found distance
1588	dist_dx = CEILING( l_wall(k,j,i) / dx )
1589
1590	ENDDO
1591
1592	ENDIF !Check for any walls within vicinity
1593
1594	ELSE !Check if (i,j,k) belongs to atmosphere
1595
1596	l_wall(k,j,i) = l_black(k)
1597
1598	ENDIF
1599
1600	ENDDO !j loop
1601	ENDDO !i loop
1602
1603	DEALLOCATE( vicinity )
1604	DEALLOCATE( vic_yz )
1605
1606	ENDIF !check vertical size of vicinity
1607
1608	ENDDO !k loop
1609
1610	!$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
1611
1612	ENDIF !LES or RANS mode
1613
1614	!
1615	!-- Set lateral boundary conditions for l_wall
1616	CALL exchange_horiz( l_wall, nbgp )
1617
1618	!$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
1619	!$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
1620
1621	CONTAINS
1622	!------------------------------------------------------------------------------!
1623	! Description:
1624	! ------------
1625	!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and
1626	!> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
1627	!> closest to the origin (0/0) of 'array'.
1628	!------------------------------------------------------------------------------!
1629	REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
1630
1631	IMPLICIT NONE
1632
1633	LOGICAL, INTENT(IN) :: orientation !< flag if array represents an array oriented upwards (true) or downwards (false)
1634
1635	INTEGER(iwp), INTENT(IN) :: pos_i !< x position of the yz-plane 'array'
1636
1637	INTEGER(iwp) :: a !< loop index
1638	INTEGER(iwp) :: b !< loop index
1639	INTEGER(iwp) :: jj !< loop index
1640
1641	INTEGER(KIND=1) :: maximum !< maximum of array along z dimension
1642
1643	INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
1644
1645	INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i
1646
1647	!
1648	!-- Get coordinate of first maximum along vertical dimension
1649	!-- at each y position of array (similar to function maxloc but more stable).
1650	DO a = 0, rad_j
1651	loc_k(a) = rad_k+1
1652	maximum = MAXVAL( array(:,a) )
1653	DO b = 0, rad_k+1
1654	IF ( array(b,a) == maximum ) THEN
1655	loc_k(a) = b
1656	EXIT
1657	ENDIF
1658	ENDDO
1659	ENDDO
1660	!
1661	!-- Set distance to the default maximum value (=search radius)
1662	shortest_distance = radius
1663	!
1664	!-- Calculate distance between position (0/0/0) and
1665	!-- position (pos_i/jj/loc(jj)) and only save the shortest distance.
1666	IF ( orientation ) THEN !if array is oriented upwards
1667	DO jj = 0, rad_j
1668	shortest_distance = &
1669	MIN( shortest_distance, &
1670	SQRT( MAX(REAL(pos_i, KIND=wp)dx-0.5_wpdx, 0.0_wp)**2 &
1671	+ MAX(REAL(jj, KIND=wp)dy-0.5_wpdy, 0.0_wp)**2 &
1672	+ MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2 &
1673	) &
1674	)
1675	ENDDO
1676	ELSE !if array is oriented downwards
1677	!> @note MAX within zw required to circumvent error at domain border
1678	!> At the domain border, if non-cyclic boundary is present, the
1679	!> index for zw could be -1, which will be errorneous (zw(-1) does
1680	!> not exist). The MAX function limits the index to be at least 0.
1681	DO jj = 0, rad_j
1682	shortest_distance = &
1683	MIN( shortest_distance, &
1684	SQRT( MAX(REAL(pos_i, KIND=wp)dx-0.5_wpdx, 0.0_wp)**2 &
1685	+ MAX(REAL(jj, KIND=wp)dy-0.5_wpdy, 0.0_wp)**2 &
1686	+ MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2 &
1687	) &
1688	)
1689	ENDDO
1690	ENDIF
1691
1692	END FUNCTION
1693
1694	!------------------------------------------------------------------------------!
1695	! Description:
1696	! ------------
1697	!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point
1698	!> (kp,jp,ip) containing the first bit of wall_flags_0 into the array
1699	!> 'vicinity'. Only copy first bit as this indicates the presence of topography.
1700	!------------------------------------------------------------------------------!
1701	SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir )
1702
1703	IMPLICIT NONE
1704
1705	INTEGER(iwp), INTENT(IN) :: il !< left loop boundary
1706	INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction
1707	INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary
1708	INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary
1709	INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction
1710	INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary
1711	INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary
1712	INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction
1713	INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary
1714
1715	INTEGER(iwp) :: i !< loop index
1716	INTEGER(iwp) :: j !< loop index
1717	INTEGER(iwp) :: k !< loop index
1718
1719	DO i = il, ir
1720	DO j = js, jn
1721	DO k = kb, kt
1722	vicinity(k,j,i) = MERGE( 0, 1, &
1723	BTEST( wall_flags_0(kp+k,jp+j,ip+i), 0 ) )
1724	ENDDO
1725	ENDDO
1726	ENDDO
1727
1728	END SUBROUTINE copy_into_vicinity
1729
1730	END SUBROUTINE tcm_init_mixing_length
1731
1732
1733	!------------------------------------------------------------------------------!
1734	! Description:
1735	! ------------
1736	!> Initialize virtual velocities used later in production_e.
1737	!------------------------------------------------------------------------------!
1738	SUBROUTINE production_e_init
1739
1740	USE arrays_3d, &
1741	ONLY: drho_air_zw, zu
1742
1743	USE control_parameters, &
1744	ONLY: constant_flux_layer
1745
1746	USE surface_layer_fluxes_mod, &
1747	ONLY: phi_m
1748
1749	USE surface_mod, &
1750	ONLY : surf_def_h, surf_lsm_h, surf_usm_h
1751
1752	IMPLICIT NONE
1753
1754	INTEGER(iwp) :: i !< grid index x-direction
1755	INTEGER(iwp) :: j !< grid index y-direction
1756	INTEGER(iwp) :: k !< grid index z-direction
1757	INTEGER(iwp) :: m !< running index surface elements
1758
1759	REAL(wp) :: km_sfc !< diffusion coefficient, used to compute virtual velocities
1760
1761	IF ( constant_flux_layer ) THEN
1762	!
1763	!-- Calculate a virtual velocity at the surface in a way that the
1764	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1765	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1766	!-- production term at k=1 (see production_e_ij).
1767	!-- The velocity gradient has to be limited in case of too small km
1768	!-- (otherwise the timestep may be significantly reduced by large
1769	!-- surface winds).
1770	!-- not available in case of non-cyclic boundary conditions.
1771	!-- Default surfaces, upward-facing
1772	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1773	!$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1774	!$ACC PRESENT(surf_def_h(0), u, v, drho_air_zw, zu)
1775	DO m = 1, surf_def_h(0)%ns
1776
1777	i = surf_def_h(0)%i(m)
1778	j = surf_def_h(0)%j(m)
1779	k = surf_def_h(0)%k(m)
1780	!
1781	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1782	!-- and km are not on the same grid. Actually, a further
1783	!-- interpolation of km onto the u/v-grid is necessary. However, the
1784	!-- effect of this error is negligible.
1785	km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / &
1786	phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
1787
1788	surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * &
1789	drho_air_zw(k-1) * &
1790	( zu(k+1) - zu(k-1) ) / &
1791	( km_sfc + 1.0E-20_wp )
1792	surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * &
1793	drho_air_zw(k-1) * &
1794	( zu(k+1) - zu(k-1) ) / &
1795	( km_sfc + 1.0E-20_wp )
1796
1797	IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > &
1798	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1799	) surf_def_h(0)%u_0(m) = u(k-1,j,i)
1800
1801	IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > &
1802	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1803	) surf_def_h(0)%v_0(m) = v(k-1,j,i)
1804
1805	ENDDO
1806	!
1807	!-- Default surfaces, downward-facing surfaces
1808	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1809	!$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1810	!$ACC PRESENT(surf_def_h(1), u, v, drho_air_zw, zu, km)
1811	DO m = 1, surf_def_h(1)%ns
1812
1813	i = surf_def_h(1)%i(m)
1814	j = surf_def_h(1)%j(m)
1815	k = surf_def_h(1)%k(m)
1816	!
1817	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1818	!-- and km are not on the same grid. Actually, a further
1819	!-- interpolation of km onto the u/v-grid is necessary. However, the
1820	!-- effect of this error is negligible.
1821	surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * &
1822	drho_air_zw(k-1) * &
1823	( zu(k+1) - zu(k-1) ) / &
1824	( km(k,j,i) + 1.0E-20_wp )
1825	surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * &
1826	drho_air_zw(k-1) * &
1827	( zu(k+1) - zu(k-1) ) / &
1828	( km(k,j,i) + 1.0E-20_wp )
1829
1830	IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > &
1831	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1832	) surf_def_h(1)%u_0(m) = u(k+1,j,i)
1833
1834	IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > &
1835	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1836	) surf_def_h(1)%v_0(m) = v(k+1,j,i)
1837
1838	ENDDO
1839	!
1840	!-- Natural surfaces, upward-facing
1841	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1842	!$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1843	!$ACC PRESENT(surf_lsm_h, u, v, drho_air_zw, zu)
1844	DO m = 1, surf_lsm_h%ns
1845
1846	i = surf_lsm_h%i(m)
1847	j = surf_lsm_h%j(m)
1848	k = surf_lsm_h%k(m)
1849	!
1850	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1851	!-- and km are not on the same grid. Actually, a further
1852	!-- interpolation of km onto the u/v-grid is necessary. However, the
1853	!-- effect of this error is negligible.
1854	km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / &
1855	phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
1856
1857	surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * &
1858	drho_air_zw(k-1) * &
1859	( zu(k+1) - zu(k-1) ) / &
1860	( km_sfc + 1.0E-20_wp )
1861	surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * &
1862	drho_air_zw(k-1) * &
1863	( zu(k+1) - zu(k-1) ) / &
1864	( km_sfc + 1.0E-20_wp )
1865
1866	IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > &
1867	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1868	) surf_lsm_h%u_0(m) = u(k-1,j,i)
1869
1870	IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > &
1871	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1872	) surf_lsm_h%v_0(m) = v(k-1,j,i)
1873
1874	ENDDO
1875	!
1876	!-- Urban surfaces, upward-facing
1877	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
1878	!$ACC PARALLEL LOOP PRIVATE(i, j, k, m, km_sfc) &
1879	!$ACC PRESENT(surf_usm_h, u, v, drho_air_zw, zu)
1880	DO m = 1, surf_usm_h%ns
1881
1882	i = surf_usm_h%i(m)
1883	j = surf_usm_h%j(m)
1884	k = surf_usm_h%k(m)
1885	!
1886	!-- Note, calculation of u_0 and v_0 is not fully accurate, as u/v
1887	!-- and km are not on the same grid. Actually, a further
1888	!-- interpolation of km onto the u/v-grid is necessary. However, the
1889	!-- effect of this error is negligible.
1890	km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / &
1891	phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
1892
1893	surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * &
1894	drho_air_zw(k-1) * &
1895	( zu(k+1) - zu(k-1) ) / &
1896	( km_sfc + 1.0E-20_wp )
1897	surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * &
1898	drho_air_zw(k-1) * &
1899	( zu(k+1) - zu(k-1) ) / &
1900	( km_sfc + 1.0E-20_wp )
1901
1902	IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > &
1903	ABS( u(k+1,j,i) - u(k-1,j,i) ) &
1904	) surf_usm_h%u_0(m) = u(k-1,j,i)
1905
1906	IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > &
1907	ABS( v(k+1,j,i) - v(k-1,j,i) ) &
1908	) surf_usm_h%v_0(m) = v(k-1,j,i)
1909
1910	ENDDO
1911
1912	ENDIF
1913
1914	END SUBROUTINE production_e_init
1915
1916
1917	!------------------------------------------------------------------------------!
1918	! Description:
1919	! ------------
1920	!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
1921	!> Vector-optimized version
1922	!------------------------------------------------------------------------------!
1923	SUBROUTINE tcm_prognostic_equations
1924
1925	USE arrays_3d, &
1926	ONLY: ddzu
1927
1928	USE control_parameters, &
1929	ONLY: f, scalar_advec, tsc
1930
1931	USE surface_mod, &
1932	ONLY : surf_def_h
1933
1934	IMPLICIT NONE
1935
1936	INTEGER(iwp) :: i !< loop index
1937	INTEGER(iwp) :: j !< loop index
1938	INTEGER(iwp) :: k !< loop index
1939	INTEGER(iwp) :: m !< loop index
1940	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
1941	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
1942
1943	REAL(wp) :: sbt !< wheighting factor for sub-time step
1944
1945	!
1946	!-- If required, compute prognostic equation for turbulent kinetic
1947	!-- energy (TKE)
1948	IF ( .NOT. constant_diffusion ) THEN
1949
1950	CALL cpu_log( log_point(16), 'tke-equation', 'start' )
1951
1952	sbt = tsc(2)
1953	IF ( .NOT. use_upstream_for_tke ) THEN
1954	IF ( scalar_advec == 'bc-scheme' ) THEN
1955
1956	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
1957	!
1958	!-- Bott-Chlond scheme always uses Euler time step. Thus:
1959	sbt = 1.0_wp
1960	ENDIF
1961	tend = 0.0_wp
1962	CALL advec_s_bc( e, 'e' )
1963
1964	ENDIF
1965	ENDIF
1966
1967	!
1968	!-- TKE-tendency terms with no communication
1969	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
1970	IF ( use_upstream_for_tke ) THEN
1971	tend = 0.0_wp
1972	CALL advec_s_up( e )
1973	ELSE
1974	!$ACC KERNELS PRESENT(tend)
1975	tend = 0.0_wp
1976	!$ACC END KERNELS
1977	IF ( timestep_scheme(1:5) == 'runge' ) THEN
1978	IF ( ws_scheme_sca ) THEN
1979	CALL advec_s_ws( e, 'e' )
1980	ELSE
1981	CALL advec_s_pw( e )
1982	ENDIF
1983	ELSE
1984	CALL advec_s_up( e )
1985	ENDIF
1986	ENDIF
1987	ENDIF
1988
1989	CALL production_e( .FALSE. )
1990
1991	IF ( .NOT. humidity ) THEN
1992	IF ( ocean_mode ) THEN
1993	CALL diffusion_e( prho, prho_reference )
1994	ELSE
1995	CALL diffusion_e( pt, pt_reference )
1996	ENDIF
1997	ELSE
1998	CALL diffusion_e( vpt, pt_reference )
1999	ENDIF
2000
2001	!
2002	!-- Additional sink term for flows through plant canopies
2003	IF ( plant_canopy ) CALL pcm_tendency( 6 )
2004
2005	CALL user_actions( 'e-tendency' )
2006
2007	!
2008	!-- Prognostic equation for TKE.
2009	!-- Eliminate negative TKE values, which can occur due to numerical
2010	!-- reasons in the course of the integration. In such cases the old TKE
2011	!-- value is reduced by 90%.
2012	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2013	!$ACC PRESENT(e, tend, te_m, wall_flags_0) &
2014	!$ACC PRESENT(tsc(3:3)) &
2015	!$ACC PRESENT(e_p)
2016	DO i = nxl, nxr
2017	DO j = nys, nyn
2018	DO k = nzb+1, nzt
2019	e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2020	tsc(3) * te_m(k,j,i) ) &
2021	) &
2022	* MERGE( 1.0_wp, 0.0_wp, &
2023	BTEST( wall_flags_0(k,j,i), 0 ) &
2024	)
2025	IF ( e_p(k,j,i) < 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
2026	ENDDO
2027	ENDDO
2028	ENDDO
2029
2030	!
2031	!-- Calculate tendencies for the next Runge-Kutta step
2032	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2033	IF ( intermediate_timestep_count == 1 ) THEN
2034	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2035	!$ACC PRESENT(tend, te_m)
2036	DO i = nxl, nxr
2037	DO j = nys, nyn
2038	DO k = nzb+1, nzt
2039	te_m(k,j,i) = tend(k,j,i)
2040	ENDDO
2041	ENDDO
2042	ENDDO
2043	ELSEIF ( intermediate_timestep_count < &
2044	intermediate_timestep_count_max ) THEN
2045	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
2046	!$ACC PRESENT(tend, te_m)
2047	DO i = nxl, nxr
2048	DO j = nys, nyn
2049	DO k = nzb+1, nzt
2050	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2051	+ 5.3125_wp * te_m(k,j,i)
2052	ENDDO
2053	ENDDO
2054	ENDDO
2055	ENDIF
2056	ENDIF
2057
2058	CALL cpu_log( log_point(16), 'tke-equation', 'stop' )
2059
2060	ENDIF ! TKE equation
2061
2062	!
2063	!-- If required, compute prognostic equation for TKE dissipation rate
2064	IF ( rans_tke_e ) THEN
2065
2066	CALL cpu_log( log_point(33), 'diss-equation', 'start' )
2067
2068	sbt = tsc(2)
2069	IF ( .NOT. use_upstream_for_tke ) THEN
2070	IF ( scalar_advec == 'bc-scheme' ) THEN
2071
2072	IF ( timestep_scheme(1:5) /= 'runge' ) THEN
2073	!
2074	!-- Bott-Chlond scheme always uses Euler time step. Thus:
2075	sbt = 1.0_wp
2076	ENDIF
2077	tend = 0.0_wp
2078	CALL advec_s_bc( diss, 'diss' )
2079
2080	ENDIF
2081	ENDIF
2082
2083	!
2084	!-- dissipation-tendency terms with no communication
2085	IF ( scalar_advec /= 'bc-scheme' .OR. use_upstream_for_tke ) THEN
2086	IF ( use_upstream_for_tke ) THEN
2087	tend = 0.0_wp
2088	CALL advec_s_up( diss )
2089	ELSE
2090	tend = 0.0_wp
2091	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2092	IF ( ws_scheme_sca ) THEN
2093	CALL advec_s_ws( diss, 'diss' )
2094	ELSE
2095	CALL advec_s_pw( diss )
2096	ENDIF
2097	ELSE
2098	CALL advec_s_up( diss )
2099	ENDIF
2100	ENDIF
2101	ENDIF
2102	!
2103	!-- Production of TKE dissipation rate
2104	CALL production_e( .TRUE. )
2105	!
2106	!-- Diffusion term of TKE dissipation rate
2107	CALL diffusion_diss
2108	!
2109	!-- Additional sink term for flows through plant canopies
2110	! IF ( plant_canopy ) CALL pcm_tendency( ? ) !> @todo not yet implemented
2111
2112	! CALL user_actions( 'diss-tendency' ) !> @todo not yet implemented
2113
2114	!
2115	!-- Prognostic equation for TKE dissipation.
2116	!-- Eliminate negative dissipation values, which can occur due to numerical
2117	!-- reasons in the course of the integration. In such cases the old
2118	!-- dissipation value is reduced by 90%.
2119	DO i = nxl, nxr
2120	DO j = nys, nyn
2121	DO k = nzb+1, nzt
2122	diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + &
2123	tsc(3) * tdiss_m(k,j,i) ) &
2124	) &
2125	* MERGE( 1.0_wp, 0.0_wp, &
2126	BTEST( wall_flags_0(k,j,i), 0 ) &
2127	)
2128	IF ( diss_p(k,j,i) < 0.0_wp ) &
2129	diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
2130	ENDDO
2131	ENDDO
2132	ENDDO
2133
2134	!
2135	!-- Calculate tendencies for the next Runge-Kutta step
2136	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2137	IF ( intermediate_timestep_count == 1 ) THEN
2138	DO i = nxl, nxr
2139	DO j = nys, nyn
2140	DO k = nzb+1, nzt
2141	tdiss_m(k,j,i) = tend(k,j,i)
2142	ENDDO
2143	ENDDO
2144	ENDDO
2145	ELSEIF ( intermediate_timestep_count < &
2146	intermediate_timestep_count_max ) THEN
2147	DO i = nxl, nxr
2148	DO j = nys, nyn
2149	DO k = nzb+1, nzt
2150	tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) &
2151	+ 5.3125_wp * tdiss_m(k,j,i)
2152	ENDDO
2153	ENDDO
2154	ENDDO
2155	ENDIF
2156	ENDIF
2157
2158	CALL cpu_log( log_point(33), 'diss-equation', 'stop' )
2159
2160	ENDIF
2161
2162	END SUBROUTINE tcm_prognostic_equations
2163
2164
2165	!------------------------------------------------------------------------------!
2166	! Description:
2167	! ------------
2168	!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
2169	!> Cache-optimized version
2170	!------------------------------------------------------------------------------!
2171	SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn )
2172
2173	USE arrays_3d, &
2174	ONLY: diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss, &
2175	flux_l_e, flux_s_diss, flux_s_e
2176
2177	USE control_parameters, &
2178	ONLY: tsc
2179
2180	IMPLICIT NONE
2181
2182	INTEGER(iwp) :: i !< loop index x direction
2183	INTEGER(iwp) :: i_omp !< first loop index of i-loop in prognostic_equations
2184	INTEGER(iwp) :: j !< loop index y direction
2185	INTEGER(iwp) :: k !< loop index z direction
2186	INTEGER(iwp) :: tn !< task number of openmp task
2187
2188	!
2189	!-- If required, compute prognostic equation for turbulent kinetic
2190	!-- energy (TKE)
2191	IF ( .NOT. constant_diffusion ) THEN
2192
2193	!
2194	!-- Tendency-terms for TKE
2195	tend(:,j,i) = 0.0_wp
2196	IF ( timestep_scheme(1:5) == 'runge' &
2197	.AND. .NOT. use_upstream_for_tke ) THEN
2198	IF ( ws_scheme_sca ) THEN
2199	CALL advec_s_ws( i, j, e, 'e', flux_s_e, diss_s_e, &
2200	flux_l_e, diss_l_e , i_omp, tn )
2201	ELSE
2202	CALL advec_s_pw( i, j, e )
2203	ENDIF
2204	ELSE
2205	CALL advec_s_up( i, j, e )
2206	ENDIF
2207
2208	CALL production_e( i, j, .FALSE. )
2209
2210	IF ( .NOT. humidity ) THEN
2211	IF ( ocean_mode ) THEN
2212	CALL diffusion_e( i, j, prho, prho_reference )
2213	ELSE
2214	CALL diffusion_e( i, j, pt, pt_reference )
2215	ENDIF
2216	ELSE
2217	CALL diffusion_e( i, j, vpt, pt_reference )
2218	ENDIF
2219
2220	!
2221	!-- Additional sink term for flows through plant canopies
2222	IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 )
2223
2224	CALL user_actions( i, j, 'e-tendency' )
2225
2226	!
2227	!-- Prognostic equation for TKE.
2228	!-- Eliminate negative TKE values, which can occur due to numerical
2229	!-- reasons in the course of the integration. In such cases the old
2230	!-- TKE value is reduced by 90%.
2231	DO k = nzb+1, nzt
2232	e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
2233	tsc(3) * te_m(k,j,i) ) &
2234	) &
2235	* MERGE( 1.0_wp, 0.0_wp, &
2236	BTEST( wall_flags_0(k,j,i), 0 ) &
2237	)
2238	IF ( e_p(k,j,i) <= 0.0_wp ) e_p(k,j,i) = 0.1_wp * e(k,j,i)
2239	ENDDO
2240
2241	!
2242	!-- Calculate tendencies for the next Runge-Kutta step
2243	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2244	IF ( intermediate_timestep_count == 1 ) THEN
2245	DO k = nzb+1, nzt
2246	te_m(k,j,i) = tend(k,j,i)
2247	ENDDO
2248	ELSEIF ( intermediate_timestep_count < &
2249	intermediate_timestep_count_max ) THEN
2250	DO k = nzb+1, nzt
2251	te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2252	5.3125_wp * te_m(k,j,i)
2253	ENDDO
2254	ENDIF
2255	ENDIF
2256
2257	ENDIF ! TKE equation
2258
2259	!
2260	!-- If required, compute prognostic equation for TKE dissipation rate
2261	IF ( rans_tke_e ) THEN
2262	!
2263	!-- Tendency-terms for dissipation
2264	tend(:,j,i) = 0.0_wp
2265	IF ( timestep_scheme(1:5) == 'runge' &
2266	.AND. .NOT. use_upstream_for_tke ) THEN
2267	IF ( ws_scheme_sca ) THEN
2268	CALL advec_s_ws( i, j, diss, 'diss', flux_s_diss, diss_s_diss, &
2269	flux_l_diss, diss_l_diss, i_omp, tn )
2270	ELSE
2271	CALL advec_s_pw( i, j, diss )
2272	ENDIF
2273	ELSE
2274	CALL advec_s_up( i, j, diss )
2275	ENDIF
2276	!
2277	!-- Production of TKE dissipation rate
2278	CALL production_e( i, j, .TRUE. )
2279	!
2280	!-- Diffusion term of TKE dissipation rate
2281	CALL diffusion_diss( i, j )
2282	!
2283	!-- Additional sink term for flows through plant canopies
2284	! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !> @todo not yet implemented
2285
2286	! CALL user_actions( i, j, 'diss-tendency' ) !> @todo not yet implemented
2287
2288	!
2289	!-- Prognostic equation for TKE dissipation
2290	!-- Eliminate negative dissipation values, which can occur due to
2291	!-- numerical reasons in the course of the integration. In such cases
2292	!-- the old dissipation value is reduced by 90%.
2293	DO k = nzb+1, nzt
2294	diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + &
2295	tsc(3) * tdiss_m(k,j,i) ) &
2296	) &
2297	* MERGE( 1.0_wp, 0.0_wp, &
2298	BTEST( wall_flags_0(k,j,i), 0 )&
2299	)
2300	ENDDO
2301
2302	!
2303	!-- Calculate tendencies for the next Runge-Kutta step
2304	IF ( timestep_scheme(1:5) == 'runge' ) THEN
2305	IF ( intermediate_timestep_count == 1 ) THEN
2306	DO k = nzb+1, nzt
2307	tdiss_m(k,j,i) = tend(k,j,i)
2308	ENDDO
2309	ELSEIF ( intermediate_timestep_count < &
2310	intermediate_timestep_count_max ) THEN
2311	DO k = nzb+1, nzt
2312	tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + &
2313	5.3125_wp * tdiss_m(k,j,i)
2314	ENDDO
2315	ENDIF
2316	ENDIF
2317
2318	ENDIF ! dissipation equation
2319
2320	END SUBROUTINE tcm_prognostic_equations_ij
2321
2322
2323	!------------------------------------------------------------------------------!
2324	! Description:
2325	! ------------
2326	!> Production terms (shear + buoyancy) of the TKE.
2327	!> Vector-optimized version
2328	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
2329	!> not considered well!
2330	!------------------------------------------------------------------------------!
2331	SUBROUTINE production_e( diss_production )
2332
2333	USE arrays_3d, &
2334	ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
2335
2336	USE control_parameters, &
2337	ONLY: cloud_droplets, constant_flux_layer, neutral, &
2338	rho_reference, use_single_reference_value, use_surface_fluxes, &
2339	use_top_fluxes
2340
2341	USE grid_variables, &
2342	ONLY: ddx, dx, ddy, dy
2343
2344	USE bulk_cloud_model_mod, &
2345	ONLY: bulk_cloud_model
2346
2347	USE surface_mod, &
2348	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
2349	surf_usm_v
2350
2351	IMPLICIT NONE
2352
2353	LOGICAL :: diss_production
2354
2355	INTEGER(iwp) :: i !< running index x-direction
2356	INTEGER(iwp) :: j !< running index y-direction
2357	INTEGER(iwp) :: k !< running index z-direction
2358	INTEGER(iwp) :: l !< running index for different surface type orientation
2359	INTEGER(iwp) :: m !< running index surface elements
2360	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
2361	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
2362	INTEGER(iwp) :: flag_nr !< number of masking flag
2363
2364	REAL(wp) :: def !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
2365	REAL(wp) :: flag !< flag to mask topography
2366	REAL(wp) :: k1 !< temporary factor
2367	REAL(wp) :: k2 !< temporary factor
2368	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
2369	REAL(wp) :: theta !< virtual potential temperature
2370	REAL(wp) :: temp !< theta * Exner-function
2371	REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation
2372	REAL(wp) :: usvs !< momentum flux u"v"
2373	REAL(wp) :: vsus !< momentum flux v"u"
2374	REAL(wp) :: wsus !< momentum flux w"u"
2375	REAL(wp) :: wsvs !< momentum flux w"v"
2376
2377	REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction
2378	REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction
2379	REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction
2380	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction
2381	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction
2382	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction
2383	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction
2384	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction
2385	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction
2386	REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary flux-array in z-direction
2387
2388
2389
2390	!
2391	!-- Calculate TKE production by shear. Calculate gradients at all grid
2392	!-- points first, gradients at surface-bounded grid points will be
2393	!-- overwritten further below.
2394	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) &
2395	!$ACC PRIVATE(surf_s, surf_e) &
2396	!$ACC PRIVATE(dudx(:), dudy(:), dudz(:), dvdx(:), dvdy(:), dvdz(:), dwdx(:), dwdy(:), dwdz(:)) &
2397	!$ACC PRESENT(e, u, v, w, diss, dd2zu, ddzw, km, wall_flags_0) &
2398	!$ACC PRESENT(tend) &
2399	!$ACC PRESENT(surf_def_h(0:1), surf_def_v(0:3)) &
2400	!$ACC PRESENT(surf_lsm_h, surf_lsm_v(0:3)) &
2401	!$ACC PRESENT(surf_usm_h, surf_usm_v(0:3))
2402	DO i = nxl, nxr
2403	DO j = nys, nyn
2404	!$ACC LOOP PRIVATE(k)
2405	DO k = nzb+1, nzt
2406
2407	dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
2408	dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
2409	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
2410	dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
2411	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
2412
2413	dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
2414	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
2415	dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
2416	dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
2417	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
2418
2419	dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
2420	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
2421	dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
2422	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
2423	dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
2424
2425	ENDDO
2426
2427
2428	flag_nr = 29
2429
2430
2431	IF ( constant_flux_layer ) THEN
2432	!
2433
2434	flag_nr = 0
2435
2436	!-- Position beneath wall
2437	!-- (2) - Will allways be executed.
2438	!-- 'bottom and wall: use u_0,v_0 and wall functions'
2439	!
2440	!-- Compute gradients at north- and south-facing surfaces.
2441	!-- First, for default surfaces, then for urban surfaces.
2442	!-- Note, so far no natural vertical surfaces implemented
2443	DO l = 0, 1
2444	surf_s = surf_def_v(l)%start_index(j,i)
2445	surf_e = surf_def_v(l)%end_index(j,i)
2446	!$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2447	DO m = surf_s, surf_e
2448	k = surf_def_v(l)%k(m)
2449	usvs = surf_def_v(l)%mom_flux_tke(0,m)
2450	wsvs = surf_def_v(l)%mom_flux_tke(1,m)
2451
2452	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
2453	* 0.5_wp * dy
2454	!
2455	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2456	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2457	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2458	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2459	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2460	ENDDO
2461	!
2462	!-- Natural surfaces
2463	surf_s = surf_lsm_v(l)%start_index(j,i)
2464	surf_e = surf_lsm_v(l)%end_index(j,i)
2465	!$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2466	DO m = surf_s, surf_e
2467	k = surf_lsm_v(l)%k(m)
2468	usvs = surf_lsm_v(l)%mom_flux_tke(0,m)
2469	wsvs = surf_lsm_v(l)%mom_flux_tke(1,m)
2470
2471	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
2472	* 0.5_wp * dy
2473	!
2474	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2475	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2476	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2477	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2478	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2479	ENDDO
2480	!
2481	!-- Urban surfaces
2482	surf_s = surf_usm_v(l)%start_index(j,i)
2483	surf_e = surf_usm_v(l)%end_index(j,i)
2484	!$ACC LOOP PRIVATE(m, k, usvs, wsvs, km_neutral, sign_dir)
2485	DO m = surf_s, surf_e
2486	k = surf_usm_v(l)%k(m)
2487	usvs = surf_usm_v(l)%mom_flux_tke(0,m)
2488	wsvs = surf_usm_v(l)%mom_flux_tke(1,m)
2489
2490	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
2491	* 0.5_wp * dy
2492	!
2493	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2494	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2495	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
2496	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
2497	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
2498	ENDDO
2499	ENDDO
2500	!
2501	!-- Compute gradients at east- and west-facing walls
2502	DO l = 2, 3
2503	surf_s = surf_def_v(l)%start_index(j,i)
2504	surf_e = surf_def_v(l)%end_index(j,i)
2505	!$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2506	DO m = surf_s, surf_e
2507	k = surf_def_v(l)%k(m)
2508	vsus = surf_def_v(l)%mom_flux_tke(0,m)
2509	wsus = surf_def_v(l)%mom_flux_tke(1,m)
2510
2511	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2512	* 0.5_wp * dx
2513	!
2514	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2515	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2516	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2517	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2518	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2519	ENDDO
2520	!
2521	!-- Natural surfaces
2522	surf_s = surf_lsm_v(l)%start_index(j,i)
2523	surf_e = surf_lsm_v(l)%end_index(j,i)
2524	!$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2525	DO m = surf_s, surf_e
2526	k = surf_lsm_v(l)%k(m)
2527	vsus = surf_lsm_v(l)%mom_flux_tke(0,m)
2528	wsus = surf_lsm_v(l)%mom_flux_tke(1,m)
2529
2530	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2531	* 0.5_wp * dx
2532	!
2533	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2534	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2535	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2536	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2537	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2538	ENDDO
2539	!
2540	!-- Urban surfaces
2541	surf_s = surf_usm_v(l)%start_index(j,i)
2542	surf_e = surf_usm_v(l)%end_index(j,i)
2543	!$ACC LOOP PRIVATE(m, k, vsus, wsus, km_neutral, sign_dir)
2544	DO m = surf_s, surf_e
2545	k = surf_usm_v(l)%k(m)
2546	vsus = surf_usm_v(l)%mom_flux_tke(0,m)
2547	wsus = surf_usm_v(l)%mom_flux_tke(1,m)
2548
2549	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
2550	* 0.5_wp * dx
2551	!
2552	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
2553	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
2554	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
2555	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
2556	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
2557	ENDDO
2558	ENDDO
2559	!
2560	!-- Compute gradients at upward-facing surfaces
2561	surf_s = surf_def_h(0)%start_index(j,i)
2562	surf_e = surf_def_h(0)%end_index(j,i)
2563	!$ACC LOOP PRIVATE(m, k)
2564	DO m = surf_s, surf_e
2565	k = surf_def_h(0)%k(m)
2566	!
2567	!-- Please note, actually, an interpolation of u_0 and v_0
2568	!-- onto the grid center would be required. However, this
2569	!-- would require several data transfers between 2D-grid and
2570	!-- wall type. The effect of this missing interpolation is
2571	!-- negligible. (See also production_e_init).
2572	dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
2573	dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
2574
2575	ENDDO
2576	!
2577	!-- Natural surfaces
2578	surf_s = surf_lsm_h%start_index(j,i)
2579	surf_e = surf_lsm_h%end_index(j,i)
2580	!$ACC LOOP PRIVATE(m, k)
2581	DO m = surf_s, surf_e
2582	k = surf_lsm_h%k(m)
2583
2584	dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
2585	dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
2586
2587	ENDDO
2588	!
2589	!-- Urban surfaces
2590	surf_s = surf_usm_h%start_index(j,i)
2591	surf_e = surf_usm_h%end_index(j,i)
2592	!$ACC LOOP PRIVATE(m, k)
2593	DO m = surf_s, surf_e
2594	k = surf_usm_h%k(m)
2595
2596	dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
2597	dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
2598
2599	ENDDO
2600	!
2601	!-- Compute gradients at downward-facing walls, only for
2602	!-- non-natural default surfaces
2603	surf_s = surf_def_h(1)%start_index(j,i)
2604	surf_e = surf_def_h(1)%end_index(j,i)
2605	!$ACC LOOP PRIVATE(m, k)
2606	DO m = surf_s, surf_e
2607	k = surf_def_h(1)%k(m)
2608
2609	dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
2610	dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
2611
2612	ENDDO
2613
2614
2615	ENDIF
2616
2617
2618	!$ACC LOOP PRIVATE(k, def, flag)
2619	DO k = nzb+1, nzt
2620
2621	def = 2.0_wp * ( dudx(k)2 + dvdy(k)2 + dwdz(k)**2 ) + &
2622	dudy(k)2 + dvdx(k)2 + dwdx(k)**2 + &
2623	dwdy(k)2 + dudz(k)2 + dvdz(k)**2 + &
2624	2.0_wp * ( dvdx(k)dudy(k) + dwdx(k)dudz(k) + &
2625	dwdy(k)*dvdz(k) )
2626
2627	IF ( def < 0.0_wp ) def = 0.0_wp
2628
2629	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
2630
2631	IF ( .NOT. diss_production ) THEN
2632
2633	!-- Compute tendency for TKE-production from shear
2634	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
2635
2636	ELSE
2637
2638	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
2639	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * &
2640	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
2641
2642	ENDIF
2643
2644	ENDDO
2645
2646
2647	ENDDO
2648	ENDDO
2649
2650	!
2651	!-- If required, calculate TKE production by buoyancy
2652	IF ( .NOT. neutral ) THEN
2653
2654	IF ( .NOT. humidity ) THEN
2655
2656	IF ( ocean_mode ) THEN
2657	!
2658	!-- So far in the ocean no special treatment of density flux
2659	!-- in the bottom and top surface layer
2660	DO i = nxl, nxr
2661	DO j = nys, nyn
2662
2663	DO k = nzb+1, nzt
2664	tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
2665	ENDDO
2666	!
2667	!-- Treatment of near-surface grid points, at up- and down-
2668	!-- ward facing surfaces
2669	IF ( use_surface_fluxes ) THEN
2670	DO l = 0, 1
2671	surf_s = surf_def_h(l)%start_index(j,i)
2672	surf_e = surf_def_h(l)%end_index(j,i)
2673	DO m = surf_s, surf_e
2674	k = surf_def_h(l)%k(m)
2675	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
2676	ENDDO
2677	ENDDO
2678	ENDIF
2679
2680	IF ( use_top_fluxes ) THEN
2681	surf_s = surf_def_h(2)%start_index(j,i)
2682	surf_e = surf_def_h(2)%end_index(j,i)
2683	DO m = surf_s, surf_e
2684	k = surf_def_h(2)%k(m)
2685	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
2686	ENDDO
2687	ENDIF
2688
2689	IF ( .NOT. diss_production ) THEN
2690
2691	!-- Compute tendency for TKE-production from shear
2692	DO k = nzb+1, nzt
2693	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2694	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2695	MERGE( rho_reference, prho(k,j,i), &
2696	use_single_reference_value ) )
2697	ENDDO
2698
2699	ELSE
2700
2701	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
2702	DO k = nzb+1, nzt
2703	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2704	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2705	MERGE( rho_reference, prho(k,j,i), &
2706	use_single_reference_value ) ) * &
2707	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
2708	c_3
2709	ENDDO
2710
2711	ENDIF
2712
2713	ENDDO
2714	ENDDO
2715
2716	ELSE ! or IF ( .NOT. ocean_mode ) THEN
2717
2718	!$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
2719	!$ACC PRIVATE(surf_s, surf_e) &
2720	!$ACC PRIVATE(tmp_flux(nzb+1:nzt)) &
2721	!$ACC PRESENT(e, diss, kh, pt, dd2zu, drho_air_zw, wall_flags_0) &
2722	!$ACC PRESENT(tend) &
2723	!$ACC PRESENT(surf_def_h(0:2)) &
2724	!$ACC PRESENT(surf_lsm_h) &
2725	!$ACC PRESENT(surf_usm_h)
2726	DO i = nxl, nxr
2727	DO j = nys, nyn
2728
2729	!$ACC LOOP PRIVATE(k)
2730	DO k = nzb+1, nzt
2731	tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
2732	ENDDO
2733
2734	IF ( use_surface_fluxes ) THEN
2735	!
2736	!-- Default surfaces, up- and downward-facing
2737	DO l = 0, 1
2738	surf_s = surf_def_h(l)%start_index(j,i)
2739	surf_e = surf_def_h(l)%end_index(j,i)
2740	!$ACC LOOP PRIVATE(m, k)
2741	DO m = surf_s, surf_e
2742	k = surf_def_h(l)%k(m)
2743	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
2744	ENDDO
2745	ENDDO
2746	!
2747	!-- Natural surfaces
2748	surf_s = surf_lsm_h%start_index(j,i)
2749	surf_e = surf_lsm_h%end_index(j,i)
2750	!$ACC LOOP PRIVATE(m, k)
2751	DO m = surf_s, surf_e
2752	k = surf_lsm_h%k(m)
2753	tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
2754	ENDDO
2755	!
2756	!-- Urban surfaces
2757	surf_s = surf_usm_h%start_index(j,i)
2758	surf_e = surf_usm_h%end_index(j,i)
2759	!$ACC LOOP PRIVATE(m, k)
2760	DO m = surf_s, surf_e
2761	k = surf_usm_h%k(m)
2762	tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
2763	ENDDO
2764	ENDIF
2765
2766	IF ( use_top_fluxes ) THEN
2767	surf_s = surf_def_h(2)%start_index(j,i)
2768	surf_e = surf_def_h(2)%end_index(j,i)
2769	!$ACC LOOP PRIVATE(m, k)
2770	DO m = surf_s, surf_e
2771	k = surf_def_h(2)%k(m)
2772	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
2773	ENDDO
2774	ENDIF
2775
2776	IF ( .NOT. diss_production ) THEN
2777
2778	!-- Compute tendency for TKE-production from shear
2779	!$ACC LOOP PRIVATE(k, flag)
2780	DO k = nzb+1, nzt
2781	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2782	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2783	MERGE( pt_reference, pt(k,j,i), &
2784	use_single_reference_value ) )
2785	ENDDO
2786
2787	ELSE
2788
2789	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
2790	DO k = nzb+1, nzt
2791	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2792	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2793	MERGE( pt_reference, pt(k,j,i), &
2794	use_single_reference_value ) ) * &
2795	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
2796	c_3
2797	ENDDO
2798
2799	ENDIF
2800
2801	ENDDO
2802	ENDDO
2803
2804	ENDIF ! from IF ( .NOT. ocean_mode )
2805
2806	ELSE ! or IF ( humidity ) THEN
2807
2808	DO i = nxl, nxr
2809	DO j = nys, nyn
2810
2811	DO k = nzb+1, nzt
2812
2813	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
2814	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2815	k2 = 0.61_wp * pt(k,j,i)
2816	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
2817	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
2818	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
2819	) * dd2zu(k)
2820	ELSE IF ( bulk_cloud_model ) THEN
2821	IF ( ql(k,j,i) == 0.0_wp ) THEN
2822	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2823	k2 = 0.61_wp * pt(k,j,i)
2824	ELSE
2825	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
2826	temp = theta * exner(k)
2827	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
2828	( q(k,j,i) - ql(k,j,i) ) * &
2829	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
2830	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
2831	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
2832	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
2833	ENDIF
2834	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
2835	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
2836	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
2837	) * dd2zu(k)
2838	ELSE IF ( cloud_droplets ) THEN
2839	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
2840	k2 = 0.61_wp * pt(k,j,i)
2841	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
2842	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
2843	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
2844	pt(k,j,i) * ( ql(k+1,j,i) - &
2845	ql(k-1,j,i) ) ) * dd2zu(k)
2846	ENDIF
2847
2848	ENDDO
2849
2850	IF ( use_surface_fluxes ) THEN
2851	!
2852	!-- Treat horizontal default surfaces
2853	DO l = 0, 1
2854	surf_s = surf_def_h(l)%start_index(j,i)
2855	surf_e = surf_def_h(l)%end_index(j,i)
2856	DO m = surf_s, surf_e
2857	k = surf_def_h(l)%k(m)
2858
2859	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
2860	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2861	k2 = 0.61_wp * pt(k,j,i)
2862	ELSE IF ( bulk_cloud_model ) THEN
2863	IF ( ql(k,j,i) == 0.0_wp ) THEN
2864	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2865	k2 = 0.61_wp * pt(k,j,i)
2866	ELSE
2867	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
2868	temp = theta * exner(k)
2869	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
2870	( q(k,j,i) - ql(k,j,i) ) * &
2871	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
2872	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
2873	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
2874	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
2875	ENDIF
2876	ELSE IF ( cloud_droplets ) THEN
2877	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
2878	k2 = 0.61_wp * pt(k,j,i)
2879	ENDIF
2880
2881	tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + &
2882	k2 * surf_def_h(l)%qsws(m) &
2883	) * drho_air_zw(k-1)
2884	ENDDO
2885	ENDDO
2886	!
2887	!-- Treat horizontal natural surfaces
2888	surf_s = surf_lsm_h%start_index(j,i)
2889	surf_e = surf_lsm_h%end_index(j,i)
2890	DO m = surf_s, surf_e
2891	k = surf_lsm_h%k(m)
2892
2893	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
2894	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2895	k2 = 0.61_wp * pt(k,j,i)
2896	ELSE IF ( bulk_cloud_model ) THEN
2897	IF ( ql(k,j,i) == 0.0_wp ) THEN
2898	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2899	k2 = 0.61_wp * pt(k,j,i)
2900	ELSE
2901	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
2902	temp = theta * exner(k)
2903	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
2904	( q(k,j,i) - ql(k,j,i) ) * &
2905	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
2906	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
2907	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
2908	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
2909	ENDIF
2910	ELSE IF ( cloud_droplets ) THEN
2911	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
2912	k2 = 0.61_wp * pt(k,j,i)
2913	ENDIF
2914
2915	tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + &
2916	k2 * surf_lsm_h%qsws(m) &
2917	) * drho_air_zw(k-1)
2918	ENDDO
2919	!
2920	!-- Treat horizontal urban surfaces
2921	surf_s = surf_usm_h%start_index(j,i)
2922	surf_e = surf_usm_h%end_index(j,i)
2923	DO m = surf_s, surf_e
2924	k = surf_usm_h%k(m)
2925
2926	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
2927	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2928	k2 = 0.61_wp * pt(k,j,i)
2929	ELSE IF ( bulk_cloud_model ) THEN
2930	IF ( ql(k,j,i) == 0.0_wp ) THEN
2931	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2932	k2 = 0.61_wp * pt(k,j,i)
2933	ELSE
2934	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
2935	temp = theta * exner(k)
2936	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
2937	( q(k,j,i) - ql(k,j,i) ) * &
2938	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
2939	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
2940	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
2941	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
2942	ENDIF
2943	ELSE IF ( cloud_droplets ) THEN
2944	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
2945	k2 = 0.61_wp * pt(k,j,i)
2946	ENDIF
2947
2948	tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + &
2949	k2 * surf_usm_h%qsws(m) &
2950	) * drho_air_zw(k-1)
2951	ENDDO
2952
2953	ENDIF ! from IF ( use_surface_fluxes ) THEN
2954
2955	IF ( use_top_fluxes ) THEN
2956
2957	surf_s = surf_def_h(2)%start_index(j,i)
2958	surf_e = surf_def_h(2)%end_index(j,i)
2959	DO m = surf_s, surf_e
2960	k = surf_def_h(2)%k(m)
2961
2962	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
2963	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2964	k2 = 0.61_wp * pt(k,j,i)
2965	ELSE IF ( bulk_cloud_model ) THEN
2966	IF ( ql(k,j,i) == 0.0_wp ) THEN
2967	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
2968	k2 = 0.61_wp * pt(k,j,i)
2969	ELSE
2970	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
2971	temp = theta * exner(k)
2972	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
2973	( q(k,j,i) - ql(k,j,i) ) * &
2974	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
2975	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
2976	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
2977	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
2978	ENDIF
2979	ELSE IF ( cloud_droplets ) THEN
2980	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
2981	k2 = 0.61_wp * pt(k,j,i)
2982	ENDIF
2983
2984	tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + &
2985	k2 * surf_def_h(2)%qsws(m) &
2986	) * drho_air_zw(k)
2987
2988	ENDDO
2989
2990	ENDIF ! from IF ( use_top_fluxes ) THEN
2991
2992	IF ( .NOT. diss_production ) THEN
2993
2994	!-- Compute tendency for TKE-production from shear
2995	DO k = nzb+1, nzt
2996	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
2997	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
2998	MERGE( vpt_reference, vpt(k,j,i), &
2999	use_single_reference_value ) )
3000	ENDDO
3001
3002	ELSE
3003
3004	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
3005	DO k = nzb+1, nzt
3006	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3007	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3008	MERGE( vpt_reference, vpt(k,j,i), &
3009	use_single_reference_value ) ) * &
3010	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3011	c_3
3012	ENDDO
3013
3014	ENDIF
3015
3016	ENDDO
3017	ENDDO
3018
3019	ENDIF
3020
3021	ENDIF
3022
3023	END SUBROUTINE production_e
3024
3025
3026	!------------------------------------------------------------------------------!
3027	! Description:
3028	! ------------
3029	!> Production terms (shear + buoyancy) of the TKE.
3030	!> Cache-optimized version
3031	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
3032	!> not considered well!
3033	!------------------------------------------------------------------------------!
3034	SUBROUTINE production_e_ij( i, j, diss_production )
3035
3036	USE arrays_3d, &
3037	ONLY: ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
3038
3039	USE control_parameters, &
3040	ONLY: cloud_droplets, constant_flux_layer, neutral, &
3041	rho_reference, use_single_reference_value, use_surface_fluxes, &
3042	use_top_fluxes
3043
3044	USE grid_variables, &
3045	ONLY: ddx, dx, ddy, dy
3046
3047	USE bulk_cloud_model_mod, &
3048	ONLY: bulk_cloud_model
3049
3050	USE surface_mod, &
3051	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, surf_usm_h, &
3052	surf_usm_v
3053
3054	IMPLICIT NONE
3055
3056	LOGICAL :: diss_production
3057
3058	INTEGER(iwp) :: i !< running index x-direction
3059	INTEGER(iwp) :: j !< running index y-direction
3060	INTEGER(iwp) :: k !< running index z-direction
3061	INTEGER(iwp) :: l !< running index for different surface type orientation
3062	INTEGER(iwp) :: m !< running index surface elements
3063	INTEGER(iwp) :: surf_e !< end index of surface elements at given i-j position
3064	INTEGER(iwp) :: surf_s !< start index of surface elements at given i-j position
3065	INTEGER(iwp) :: flag_nr !< number of masking flag
3066
3067	REAL(wp) :: def !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
3068	REAL(wp) :: flag !< flag to mask topography
3069	REAL(wp) :: k1 !< temporary factor
3070	REAL(wp) :: k2 !< temporary factor
3071	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
3072	REAL(wp) :: theta !< virtual potential temperature
3073	REAL(wp) :: temp !< theta * Exner-function
3074	REAL(wp) :: sign_dir !< sign of wall-tke flux, depending on wall orientation
3075	REAL(wp) :: usvs !< momentum flux u"v"
3076	REAL(wp) :: vsus !< momentum flux v"u"
3077	REAL(wp) :: wsus !< momentum flux w"u"
3078	REAL(wp) :: wsvs !< momentum flux w"v"
3079
3080	REAL(wp), DIMENSION(nzb+1:nzt) :: dudx !< Gradient of u-component in x-direction
3081	REAL(wp), DIMENSION(nzb+1:nzt) :: dudy !< Gradient of u-component in y-direction
3082	REAL(wp), DIMENSION(nzb+1:nzt) :: dudz !< Gradient of u-component in z-direction
3083	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdx !< Gradient of v-component in x-direction
3084	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdy !< Gradient of v-component in y-direction
3085	REAL(wp), DIMENSION(nzb+1:nzt) :: dvdz !< Gradient of v-component in z-direction
3086	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdx !< Gradient of w-component in x-direction
3087	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdy !< Gradient of w-component in y-direction
3088	REAL(wp), DIMENSION(nzb+1:nzt) :: dwdz !< Gradient of w-component in z-direction
3089	REAL(wp), DIMENSION(nzb+1:nzt) :: tmp_flux !< temporary flux-array in z-direction
3090
3091
3092
3093	!
3094	!-- Calculate TKE production by shear. Calculate gradients at all grid
3095	!-- points first, gradients at surface-bounded grid points will be
3096	!-- overwritten further below.
3097	DO k = nzb+1, nzt
3098
3099	dudx(k) = ( u(k,j,i+1) - u(k,j,i) ) * ddx
3100	dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
3101	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
3102	dudz(k) = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
3103	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
3104
3105	dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
3106	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
3107	dvdy(k) = ( v(k,j+1,i) - v(k,j,i) ) * ddy
3108	dvdz(k) = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
3109	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
3110
3111	dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
3112	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
3113	dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
3114	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
3115	dwdz(k) = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
3116
3117	ENDDO
3118
3119	flag_nr = 29
3120
3121	IF ( constant_flux_layer ) THEN
3122
3123	flag_nr = 0
3124
3125	!-- Position beneath wall
3126	!-- (2) - Will allways be executed.
3127	!-- 'bottom and wall: use u_0,v_0 and wall functions'
3128	!
3129	!-- Compute gradients at north- and south-facing surfaces.
3130	!-- First, for default surfaces, then for urban surfaces.
3131	!-- Note, so far no natural vertical surfaces implemented
3132	DO l = 0, 1
3133	surf_s = surf_def_v(l)%start_index(j,i)
3134	surf_e = surf_def_v(l)%end_index(j,i)
3135	DO m = surf_s, surf_e
3136	k = surf_def_v(l)%k(m)
3137	usvs = surf_def_v(l)%mom_flux_tke(0,m)
3138	wsvs = surf_def_v(l)%mom_flux_tke(1,m)
3139
3140	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3141	* 0.5_wp * dy
3142	!
3143	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3144	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3145	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3146	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3147	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3148	ENDDO
3149	!
3150	!-- Natural surfaces
3151	surf_s = surf_lsm_v(l)%start_index(j,i)
3152	surf_e = surf_lsm_v(l)%end_index(j,i)
3153	DO m = surf_s, surf_e
3154	k = surf_lsm_v(l)%k(m)
3155	usvs = surf_lsm_v(l)%mom_flux_tke(0,m)
3156	wsvs = surf_lsm_v(l)%mom_flux_tke(1,m)
3157
3158	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3159	* 0.5_wp * dy
3160	!
3161	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3162	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3163	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3164	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3165	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3166	ENDDO
3167	!
3168	!-- Urban surfaces
3169	surf_s = surf_usm_v(l)%start_index(j,i)
3170	surf_e = surf_usm_v(l)%end_index(j,i)
3171	DO m = surf_s, surf_e
3172	k = surf_usm_v(l)%k(m)
3173	usvs = surf_usm_v(l)%mom_flux_tke(0,m)
3174	wsvs = surf_usm_v(l)%mom_flux_tke(1,m)
3175
3176	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
3177	* 0.5_wp * dy
3178	!
3179	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3180	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3181	BTEST( wall_flags_0(k,j-1,i), flag_nr ) )
3182	dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
3183	dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
3184	ENDDO
3185	ENDDO
3186	!
3187	!-- Compute gradients at east- and west-facing walls
3188	DO l = 2, 3
3189	surf_s = surf_def_v(l)%start_index(j,i)
3190	surf_e = surf_def_v(l)%end_index(j,i)
3191	DO m = surf_s, surf_e
3192	k = surf_def_v(l)%k(m)
3193	vsus = surf_def_v(l)%mom_flux_tke(0,m)
3194	wsus = surf_def_v(l)%mom_flux_tke(1,m)
3195
3196	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3197	* 0.5_wp * dx
3198	!
3199	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3200	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3201	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3202	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3203	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3204	ENDDO
3205	!
3206	!-- Natural surfaces
3207	surf_s = surf_lsm_v(l)%start_index(j,i)
3208	surf_e = surf_lsm_v(l)%end_index(j,i)
3209	DO m = surf_s, surf_e
3210	k = surf_lsm_v(l)%k(m)
3211	vsus = surf_lsm_v(l)%mom_flux_tke(0,m)
3212	wsus = surf_lsm_v(l)%mom_flux_tke(1,m)
3213
3214	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3215	* 0.5_wp * dx
3216	!
3217	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3218	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3219	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3220	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3221	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3222	ENDDO
3223	!
3224	!-- Urban surfaces
3225	surf_s = surf_usm_v(l)%start_index(j,i)
3226	surf_e = surf_usm_v(l)%end_index(j,i)
3227	DO m = surf_s, surf_e
3228	k = surf_usm_v(l)%k(m)
3229	vsus = surf_usm_v(l)%mom_flux_tke(0,m)
3230	wsus = surf_usm_v(l)%mom_flux_tke(1,m)
3231
3232	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
3233	* 0.5_wp * dx
3234	!
3235	!-- -1.0 for right-facing wall, 1.0 for left-facing wall
3236	sign_dir = MERGE( 1.0_wp, -1.0_wp, &
3237	BTEST( wall_flags_0(k,j,i-1), flag_nr ) )
3238	dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
3239	dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
3240	ENDDO
3241	ENDDO
3242	!
3243	!-- Compute gradients at upward-facing surfaces
3244	surf_s = surf_def_h(0)%start_index(j,i)
3245	surf_e = surf_def_h(0)%end_index(j,i)
3246	DO m = surf_s, surf_e
3247	k = surf_def_h(0)%k(m)
3248	!
3249	!-- Please note, actually, an interpolation of u_0 and v_0
3250	!-- onto the grid center would be required. However, this
3251	!-- would require several data transfers between 2D-grid and
3252	!-- wall type. The effect of this missing interpolation is
3253	!-- negligible. (See also production_e_init).
3254	dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
3255	dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
3256
3257	ENDDO
3258	!
3259	!-- Natural surfaces
3260	surf_s = surf_lsm_h%start_index(j,i)
3261	surf_e = surf_lsm_h%end_index(j,i)
3262	DO m = surf_s, surf_e
3263	k = surf_lsm_h%k(m)
3264
3265	dudz(k) = ( u(k+1,j,i) - surf_lsm_h%u_0(m) ) * dd2zu(k)
3266	dvdz(k) = ( v(k+1,j,i) - surf_lsm_h%v_0(m) ) * dd2zu(k)
3267
3268	ENDDO
3269	!
3270	!-- Urban surfaces
3271	surf_s = surf_usm_h%start_index(j,i)
3272	surf_e = surf_usm_h%end_index(j,i)
3273	DO m = surf_s, surf_e
3274	k = surf_usm_h%k(m)
3275
3276	dudz(k) = ( u(k+1,j,i) - surf_usm_h%u_0(m) ) * dd2zu(k)
3277	dvdz(k) = ( v(k+1,j,i) - surf_usm_h%v_0(m) ) * dd2zu(k)
3278
3279	ENDDO
3280	!
3281	!-- Compute gradients at downward-facing walls, only for
3282	!-- non-natural default surfaces
3283	surf_s = surf_def_h(1)%start_index(j,i)
3284	surf_e = surf_def_h(1)%end_index(j,i)
3285	DO m = surf_s, surf_e
3286	k = surf_def_h(1)%k(m)
3287
3288	dudz(k) = ( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) * dd2zu(k)
3289	dvdz(k) = ( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) * dd2zu(k)
3290
3291	ENDDO
3292
3293	ENDIF
3294
3295	DO k = nzb+1, nzt
3296
3297	def = 2.0_wp * ( dudx(k)2 + dvdy(k)2 + dwdz(k)**2 ) + &
3298	dudy(k)2 + dvdx(k)2 + dwdx(k)**2 + &
3299	dwdy(k)2 + dudz(k)2 + dvdz(k)**2 + &
3300	2.0_wp * ( dvdx(k)dudy(k) + dwdx(k)dudz(k) + &
3301	dwdy(k)*dvdz(k) )
3302
3303	IF ( def < 0.0_wp ) def = 0.0_wp
3304
3305	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),flag_nr) )
3306
3307	IF ( .NOT. diss_production ) THEN
3308
3309	!-- Compute tendency for TKE-production from shear
3310	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
3311
3312	ELSE
3313
3314	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
3315	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * &
3316	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
3317
3318	ENDIF
3319
3320	ENDDO
3321
3322	!
3323	!-- If required, calculate TKE production by buoyancy
3324	IF ( .NOT. neutral ) THEN
3325
3326	IF ( .NOT. humidity ) THEN
3327
3328	IF ( ocean_mode ) THEN
3329	!
3330	!-- So far in the ocean no special treatment of density flux
3331	!-- in the bottom and top surface layer
3332	DO k = nzb+1, nzt
3333	tmp_flux(k) = kh(k,j,i) * ( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
3334	ENDDO
3335	!
3336	!-- Treatment of near-surface grid points, at up- and down-
3337	!-- ward facing surfaces
3338	IF ( use_surface_fluxes ) THEN
3339	DO l = 0, 1
3340	surf_s = surf_def_h(l)%start_index(j,i)
3341	surf_e = surf_def_h(l)%end_index(j,i)
3342	DO m = surf_s, surf_e
3343	k = surf_def_h(l)%k(m)
3344	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3345	ENDDO
3346	ENDDO
3347	ENDIF
3348
3349	IF ( use_top_fluxes ) THEN
3350	surf_s = surf_def_h(2)%start_index(j,i)
3351	surf_e = surf_def_h(2)%end_index(j,i)
3352	DO m = surf_s, surf_e
3353	k = surf_def_h(2)%k(m)
3354	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3355	ENDDO
3356	ENDIF
3357
3358	IF ( .NOT. diss_production ) THEN
3359
3360	!-- Compute tendency for TKE-production from shear
3361	DO k = nzb+1, nzt
3362	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3363	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3364	MERGE( rho_reference, prho(k,j,i), &
3365	use_single_reference_value ) )
3366	ENDDO
3367
3368	ELSE
3369
3370	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
3371	DO k = nzb+1, nzt
3372	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3373	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3374	MERGE( rho_reference, prho(k,j,i), &
3375	use_single_reference_value ) ) * &
3376	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3377	c_3
3378	ENDDO
3379
3380	ENDIF
3381
3382
3383	ELSE ! or IF ( .NOT. ocean_mode ) THEN
3384
3385	DO k = nzb+1, nzt
3386	tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
3387	ENDDO
3388
3389	IF ( use_surface_fluxes ) THEN
3390	!
3391	!-- Default surfaces, up- and downward-facing
3392	DO l = 0, 1
3393	surf_s = surf_def_h(l)%start_index(j,i)
3394	surf_e = surf_def_h(l)%end_index(j,i)
3395	DO m = surf_s, surf_e
3396	k = surf_def_h(l)%k(m)
3397	tmp_flux(k) = drho_air_zw(k-1) * surf_def_h(l)%shf(m)
3398	ENDDO
3399	ENDDO
3400	!
3401	!-- Natural surfaces
3402	surf_s = surf_lsm_h%start_index(j,i)
3403	surf_e = surf_lsm_h%end_index(j,i)
3404	DO m = surf_s, surf_e
3405	k = surf_lsm_h%k(m)
3406	tmp_flux(k) = drho_air_zw(k-1) * surf_lsm_h%shf(m)
3407	ENDDO
3408	!
3409	!-- Urban surfaces
3410	surf_s = surf_usm_h%start_index(j,i)
3411	surf_e = surf_usm_h%end_index(j,i)
3412	DO m = surf_s, surf_e
3413	k = surf_usm_h%k(m)
3414	tmp_flux(k) = drho_air_zw(k-1) * surf_usm_h%shf(m)
3415	ENDDO
3416	ENDIF
3417
3418	IF ( use_top_fluxes ) THEN
3419	surf_s = surf_def_h(2)%start_index(j,i)
3420	surf_e = surf_def_h(2)%end_index(j,i)
3421	DO m = surf_s, surf_e
3422	k = surf_def_h(2)%k(m)
3423	tmp_flux(k) = drho_air_zw(k) * surf_def_h(2)%shf(m)
3424	ENDDO
3425	ENDIF
3426
3427	IF ( .NOT. diss_production ) THEN
3428
3429	!-- Compute tendency for TKE-production from shear
3430	DO k = nzb+1, nzt
3431	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3432	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3433	MERGE( pt_reference, pt(k,j,i), &
3434	use_single_reference_value ) )
3435	ENDDO
3436
3437	ELSE
3438
3439	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
3440	DO k = nzb+1, nzt
3441	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3442	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3443	MERGE( pt_reference, pt(k,j,i), &
3444	use_single_reference_value ) ) * &
3445	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3446	c_3
3447	ENDDO
3448
3449	ENDIF
3450
3451	ENDIF ! from IF ( .NOT. ocean_mode )
3452
3453	ELSE ! or IF ( humidity ) THEN
3454
3455	DO k = nzb+1, nzt
3456
3457	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3458	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3459	k2 = 0.61_wp * pt(k,j,i)
3460	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3461	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3462	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
3463	) * dd2zu(k)
3464	ELSE IF ( bulk_cloud_model ) THEN
3465	IF ( ql(k,j,i) == 0.0_wp ) THEN
3466	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3467	k2 = 0.61_wp * pt(k,j,i)
3468	ELSE
3469	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3470	temp = theta * exner(k)
3471	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3472	( q(k,j,i) - ql(k,j,i) ) * &
3473	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3474	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3475	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3476	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3477	ENDIF
3478	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3479	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3480	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
3481	) * dd2zu(k)
3482	ELSE IF ( cloud_droplets ) THEN
3483	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3484	k2 = 0.61_wp * pt(k,j,i)
3485	tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
3486	( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) + &
3487	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
3488	pt(k,j,i) * ( ql(k+1,j,i) - &
3489	ql(k-1,j,i) ) ) * dd2zu(k)
3490	ENDIF
3491
3492	ENDDO
3493
3494	IF ( use_surface_fluxes ) THEN
3495	!
3496	!-- Treat horizontal default surfaces
3497	DO l = 0, 1
3498	surf_s = surf_def_h(l)%start_index(j,i)
3499	surf_e = surf_def_h(l)%end_index(j,i)
3500	DO m = surf_s, surf_e
3501	k = surf_def_h(l)%k(m)
3502
3503	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3504	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3505	k2 = 0.61_wp * pt(k,j,i)
3506	ELSE IF ( bulk_cloud_model ) THEN
3507	IF ( ql(k,j,i) == 0.0_wp ) THEN
3508	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3509	k2 = 0.61_wp * pt(k,j,i)
3510	ELSE
3511	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3512	temp = theta * exner(k)
3513	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3514	( q(k,j,i) - ql(k,j,i) ) * &
3515	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3516	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3517	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3518	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3519	ENDIF
3520	ELSE IF ( cloud_droplets ) THEN
3521	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3522	k2 = 0.61_wp * pt(k,j,i)
3523	ENDIF
3524
3525	tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + &
3526	k2 * surf_def_h(l)%qsws(m) &
3527	) * drho_air_zw(k-1)
3528	ENDDO
3529	ENDDO
3530	!
3531	!-- Treat horizontal natural surfaces
3532	surf_s = surf_lsm_h%start_index(j,i)
3533	surf_e = surf_lsm_h%end_index(j,i)
3534	DO m = surf_s, surf_e
3535	k = surf_lsm_h%k(m)
3536
3537	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3538	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3539	k2 = 0.61_wp * pt(k,j,i)
3540	ELSE IF ( bulk_cloud_model ) THEN
3541	IF ( ql(k,j,i) == 0.0_wp ) THEN
3542	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3543	k2 = 0.61_wp * pt(k,j,i)
3544	ELSE
3545	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3546	temp = theta * exner(k)
3547	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3548	( q(k,j,i) - ql(k,j,i) ) * &
3549	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3550	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3551	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3552	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3553	ENDIF
3554	ELSE IF ( cloud_droplets ) THEN
3555	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3556	k2 = 0.61_wp * pt(k,j,i)
3557	ENDIF
3558
3559	tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + &
3560	k2 * surf_lsm_h%qsws(m) &
3561	) * drho_air_zw(k-1)
3562	ENDDO
3563	!
3564	!-- Treat horizontal urban surfaces
3565	surf_s = surf_usm_h%start_index(j,i)
3566	surf_e = surf_usm_h%end_index(j,i)
3567	DO m = surf_s, surf_e
3568	k = surf_usm_h%k(m)
3569
3570	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3571	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3572	k2 = 0.61_wp * pt(k,j,i)
3573	ELSE IF ( bulk_cloud_model ) THEN
3574	IF ( ql(k,j,i) == 0.0_wp ) THEN
3575	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3576	k2 = 0.61_wp * pt(k,j,i)
3577	ELSE
3578	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3579	temp = theta * exner(k)
3580	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3581	( q(k,j,i) - ql(k,j,i) ) * &
3582	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3583	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3584	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3585	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3586	ENDIF
3587	ELSE IF ( cloud_droplets ) THEN
3588	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3589	k2 = 0.61_wp * pt(k,j,i)
3590	ENDIF
3591
3592	tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + &
3593	k2 * surf_usm_h%qsws(m) &
3594	) * drho_air_zw(k-1)
3595	ENDDO
3596
3597	ENDIF ! from IF ( use_surface_fluxes ) THEN
3598
3599	IF ( use_top_fluxes ) THEN
3600
3601	surf_s = surf_def_h(2)%start_index(j,i)
3602	surf_e = surf_def_h(2)%end_index(j,i)
3603	DO m = surf_s, surf_e
3604	k = surf_def_h(2)%k(m)
3605
3606	IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
3607	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3608	k2 = 0.61_wp * pt(k,j,i)
3609	ELSE IF ( bulk_cloud_model ) THEN
3610	IF ( ql(k,j,i) == 0.0_wp ) THEN
3611	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
3612	k2 = 0.61_wp * pt(k,j,i)
3613	ELSE
3614	theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
3615	temp = theta * exner(k)
3616	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
3617	( q(k,j,i) - ql(k,j,i) ) * &
3618	( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) / &
3619	( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp * &
3620	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
3621	k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
3622	ENDIF
3623	ELSE IF ( cloud_droplets ) THEN
3624	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
3625	k2 = 0.61_wp * pt(k,j,i)
3626	ENDIF
3627
3628	tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + &
3629	k2 * surf_def_h(2)%qsws(m) &
3630	) * drho_air_zw(k)
3631
3632	ENDDO
3633
3634	ENDIF ! from IF ( use_top_fluxes ) THEN
3635
3636	IF ( .NOT. diss_production ) THEN
3637
3638	!-- Compute tendency for TKE-production from shear
3639	DO k = nzb+1, nzt
3640	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3641	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3642	MERGE( vpt_reference, vpt(k,j,i), &
3643	use_single_reference_value ) )
3644	ENDDO
3645
3646	ELSE
3647
3648	!-- RANS mode: Compute tendency for dissipation-rate-production from shear
3649	DO k = nzb+1, nzt
3650	flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_0(k,j,i),0) )
3651	tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
3652	MERGE( vpt_reference, vpt(k,j,i), &
3653	use_single_reference_value ) ) * &
3654	diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * &
3655	c_3
3656	ENDDO
3657
3658	ENDIF
3659
3660	ENDIF
3661
3662	ENDIF
3663
3664	END SUBROUTINE production_e_ij
3665
3666
3667	!------------------------------------------------------------------------------!
3668	! Description:
3669	! ------------
3670	!> Diffusion and dissipation terms for the TKE.
3671	!> Vector-optimized version
3672	!------------------------------------------------------------------------------!
3673	SUBROUTINE diffusion_e( var, var_reference )
3674
3675	USE arrays_3d, &
3676	ONLY: ddzu, dd2zu, ddzw, drho_air, rho_air_zw
3677
3678	USE control_parameters, &
3679	ONLY: atmos_ocean_sign, use_single_reference_value, &
3680	wall_adjustment, wall_adjustment_factor
3681
3682	USE grid_variables, &
3683	ONLY: ddx2, ddy2
3684
3685	USE bulk_cloud_model_mod, &
3686	ONLY: collision_turbulence
3687
3688	USE particle_attributes, &
3689	ONLY: use_sgs_for_particles, wang_kernel
3690
3691	USE surface_mod, &
3692	ONLY : bc_h
3693
3694	IMPLICIT NONE
3695
3696	INTEGER(iwp) :: i !< running index x direction
3697	INTEGER(iwp) :: j !< running index y direction
3698	INTEGER(iwp) :: k !< running index z direction
3699	INTEGER(iwp) :: m !< running index surface elements
3700
3701	REAL(wp) :: flag !< flag to mask topography
3702	REAL(wp) :: l !< mixing length
3703	REAL(wp) :: ll !< adjusted l
3704	REAL(wp) :: var_reference !< reference temperature
3705	#ifdef _OPENACC
3706	!
3707	!-- From mixing_length_les:
3708	REAL(wp) :: l_stable !< mixing length according to stratification
3709	!
3710	!-- From mixing_length_rans:
3711	REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector
3712	REAL(wp) :: l_diss !< mixing length for dissipation
3713	REAL(wp) :: rif !< Richardson flux number
3714	!
3715	!-- From both:
3716	REAL(wp) :: dvar_dz !< vertical gradient of var
3717	#endif
3718
3719	REAL(wp) :: dissipation !< TKE dissipation
3720
3721	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
3722
3723
3724	!
3725	!-- Calculate the tendency terms
3726	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
3727	!$ACC PRIVATE(flag, l, ll, dissipation) &
3728	!$ACC PRIVATE(l_stable, duv2_dz2, l_diss, rif, dvar_dz) &
3729	!$ACC PRESENT(e, u, v, km, var, wall_flags_0) &
3730	!$ACC PRESENT(ddzu, dd2zu, ddzw, rho_air_zw, drho_air) &
3731	!$ACC PRESENT(l_black, l_grid, l_wall) &
3732	!$ACC PRESENT(diss, tend)
3733	DO i = nxl, nxr
3734	DO j = nys, nyn
3735	DO k = nzb+1, nzt
3736	!
3737	!-- Predetermine flag to mask topography
3738	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3739
3740	!
3741	!-- Calculate dissipation
3742	IF ( les_dynamic .OR. les_mw ) THEN
3743
3744	#ifdef _OPENACC
3745	!
3746	!-- Cannot call subroutine mixing_length_les because after adding all required
3747	!-- OpenACC directives the execution crashed reliably when inside the called
3748	!-- subroutine. I'm not sure why that is...
3749	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
3750	IF ( dvar_dz > 0.0_wp ) THEN
3751	IF ( use_single_reference_value ) THEN
3752	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
3753	/ SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
3754	ELSE
3755	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
3756	/ SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
3757	ENDIF
3758	ELSE
3759	l_stable = l_grid(k)
3760	ENDIF
3761	!
3762	!-- Adjustment of the mixing length
3763	IF ( wall_adjustment ) THEN
3764	l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
3765	ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
3766	ELSE
3767	l = MIN( l_grid(k), l_stable )
3768	ll = l_grid(k)
3769	ENDIF
3770	#else
3771	CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
3772	#endif
3773
3774	dissipation = ( 0.19_wp + 0.74_wp * l / ll ) &
3775	* e(k,j,i) * SQRT( e(k,j,i) ) / l
3776
3777	ELSEIF ( rans_tke_l ) THEN
3778
3779	#ifdef _OPENACC
3780	!
3781	!-- Same reason as above...
3782	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
3783
3784	duv2_dz2 = ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2 &
3785	+ ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 &
3786	+ 1E-30_wp
3787
3788	IF ( use_single_reference_value ) THEN
3789	rif = g / var_reference * dvar_dz / duv2_dz2
3790	ELSE
3791	rif = g / var(k,j,i) * dvar_dz / duv2_dz2
3792	ENDIF
3793
3794	rif = MAX( rif, -5.0_wp )
3795	rif = MIN( rif, 1.0_wp )
3796
3797	!
3798	!-- Calculate diabatic mixing length using Dyer-profile functions
3799	IF ( rif >= 0.0_wp ) THEN
3800	l = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
3801	l_diss = l
3802	ELSE
3803	!
3804	!-- In case of unstable stratification, use mixing length of neutral case
3805	!-- for l, but consider profile functions for l_diss
3806	l = l_wall(k,j,i)
3807	l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
3808	ENDIF
3809	#else
3810	CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
3811	#endif
3812
3813	dissipation = c_0*3 e(k,j,i) * SQRT( e(k,j,i) ) / ll
3814
3815	diss(k,j,i) = dissipation * flag
3816
3817	ELSEIF ( rans_tke_e ) THEN
3818
3819	dissipation = diss(k,j,i)
3820
3821	ENDIF
3822
3823	tend(k,j,i) = tend(k,j,i) + ( &
3824	( &
3825	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
3826	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
3827	) * ddx2 * flag &
3828	+ ( &
3829	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
3830	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
3831	) * ddy2 * flag &
3832	+ ( &
3833	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
3834	* rho_air_zw(k) &
3835	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
3836	* rho_air_zw(k-1) &
3837	) * ddzw(k) * drho_air(k) &
3838	) * flag * dsig_e &
3839	- dissipation * flag
3840	!
3841	!-- Store dissipation if needed for calculating the sgs particle
3842	!-- velocities
3843	IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. &
3844	wang_kernel .OR. collision_turbulence ) ) THEN
3845	diss(k,j,i) = dissipation * flag
3846	ENDIF
3847
3848	ENDDO
3849	ENDDO
3850	ENDDO
3851
3852	!
3853	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
3854	IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. &
3855	wang_kernel .OR. collision_turbulence ) ) THEN
3856	!
3857	!-- Upward facing surfaces
3858	DO m = 1, bc_h(0)%ns
3859	i = bc_h(0)%i(m)
3860	j = bc_h(0)%j(m)
3861	k = bc_h(0)%k(m)
3862	diss(k-1,j,i) = diss(k,j,i)
3863	ENDDO
3864	!
3865	!-- Downward facing surfaces
3866	DO m = 1, bc_h(1)%ns
3867	i = bc_h(1)%i(m)
3868	j = bc_h(1)%j(m)
3869	k = bc_h(1)%k(m)
3870	diss(k+1,j,i) = diss(k,j,i)
3871	ENDDO
3872
3873	ENDIF
3874
3875	END SUBROUTINE diffusion_e
3876
3877
3878	!------------------------------------------------------------------------------!
3879	! Description:
3880	! ------------
3881	!> Diffusion and dissipation terms for the TKE.
3882	!> Cache-optimized version
3883	!------------------------------------------------------------------------------!
3884	SUBROUTINE diffusion_e_ij( i, j, var, var_reference )
3885
3886	USE arrays_3d, &
3887	ONLY: ddzu, ddzw, drho_air, rho_air_zw
3888
3889	USE grid_variables, &
3890	ONLY: ddx2, ddy2
3891
3892	USE bulk_cloud_model_mod, &
3893	ONLY: collision_turbulence
3894
3895	USE particle_attributes, &
3896	ONLY: use_sgs_for_particles, wang_kernel
3897
3898	USE surface_mod, &
3899	ONLY : bc_h
3900
3901	IMPLICIT NONE
3902
3903	INTEGER(iwp) :: i !< running index x direction
3904	INTEGER(iwp) :: j !< running index y direction
3905	INTEGER(iwp) :: k !< running index z direction
3906	INTEGER(iwp) :: m !< running index surface elements
3907	INTEGER(iwp) :: surf_e !< End index of surface elements at (j,i)-gridpoint
3908	INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint
3909
3910	REAL(wp) :: flag !< flag to mask topography
3911	REAL(wp) :: l !< mixing length
3912	REAL(wp) :: ll !< adjusted l
3913	REAL(wp) :: var_reference !< reference temperature
3914
3915	REAL(wp), DIMENSION(nzb+1:nzt) :: dissipation !< dissipation of TKE
3916
3917	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
3918
3919	!
3920	!-- Calculate the mixing length (for dissipation)
3921	DO k = nzb+1, nzt
3922	!
3923	!-- Predetermine flag to mask topography
3924	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
3925
3926	!
3927	!-- Calculate dissipation...
3928	!-- ...in case of LES
3929	IF ( les_dynamic .OR. les_mw ) THEN
3930
3931	CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
3932
3933	dissipation(k) = ( 0.19_wp + 0.74_wp * l / ll ) &
3934	* e(k,j,i) * SQRT( e(k,j,i) ) / l
3935
3936	!
3937	!-- ...in case of RANS
3938	ELSEIF ( rans_tke_l ) THEN
3939
3940	CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
3941
3942	dissipation(k) = c_0*3 e(k,j,i) * SQRT( e(k,j,i) ) / ll
3943
3944	diss(k,j,i) = dissipation(k) * flag
3945
3946	ELSEIF ( rans_tke_e ) THEN
3947
3948	dissipation(k) = diss(k,j,i)
3949
3950	ENDIF
3951
3952	!
3953	!-- Calculate the tendency term
3954	tend(k,j,i) = tend(k,j,i) + ( &
3955	( &
3956	( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) ) &
3957	- ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) ) &
3958	) * ddx2 &
3959	+ ( &
3960	( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) ) &
3961	- ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) ) &
3962	) * ddy2 &
3963	+ ( &
3964	( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1) &
3965	* rho_air_zw(k) &
3966	- ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k) &
3967	* rho_air_zw(k-1) &
3968	) * ddzw(k) * drho_air(k) &
3969	) * flag * dsig_e &
3970	- dissipation(k) * flag
3971
3972	ENDDO
3973
3974	!
3975	!-- Store dissipation if needed for calculating the sgs particle velocities
3976	IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. wang_kernel &
3977	.OR. collision_turbulence ) ) THEN
3978	DO k = nzb+1, nzt
3979	diss(k,j,i) = dissipation(k) * MERGE( 1.0_wp, 0.0_wp, &
3980	BTEST( wall_flags_0(k,j,i), 0 ) )
3981	ENDDO
3982	!
3983	!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
3984	!-- For each surface type determine start and end index (in case of elevated
3985	!-- topography several up/downward facing surfaces may exist.
3986	surf_s = bc_h(0)%start_index(j,i)
3987	surf_e = bc_h(0)%end_index(j,i)
3988	DO m = surf_s, surf_e
3989	k = bc_h(0)%k(m)
3990	diss(k-1,j,i) = diss(k,j,i)
3991	ENDDO
3992	!
3993	!-- Downward facing surfaces
3994	surf_s = bc_h(1)%start_index(j,i)
3995	surf_e = bc_h(1)%end_index(j,i)
3996	DO m = surf_s, surf_e
3997	k = bc_h(1)%k(m)
3998	diss(k+1,j,i) = diss(k,j,i)
3999	ENDDO
4000	ENDIF
4001
4002	END SUBROUTINE diffusion_e_ij
4003
4004
4005	!------------------------------------------------------------------------------!
4006	! Description:
4007	! ------------
4008	!> Diffusion term for the TKE dissipation rate
4009	!> Vector-optimized version
4010	!------------------------------------------------------------------------------!
4011	SUBROUTINE diffusion_diss()
4012	USE arrays_3d, &
4013	ONLY: ddzu, ddzw, drho_air, rho_air_zw
4014
4015	USE grid_variables, &
4016	ONLY: ddx2, ddy2
4017
4018	IMPLICIT NONE
4019
4020	INTEGER(iwp) :: i !< running index x direction
4021	INTEGER(iwp) :: j !< running index y direction
4022	INTEGER(iwp) :: k !< running index z direction
4023
4024	REAL(wp) :: flag !< flag to mask topography
4025
4026	!
4027	!-- Calculate the tendency terms
4028	DO i = nxl, nxr
4029	DO j = nys, nyn
4030	DO k = nzb+1, nzt
4031
4032	!
4033	!-- Predetermine flag to mask topography
4034	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4035
4036	tend(k,j,i) = tend(k,j,i) + &
4037	( ( &
4038	( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) &
4039	- ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) &
4040	) * ddx2 &
4041	+ ( &
4042	( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) &
4043	- ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) &
4044	) * ddy2 &
4045	+ ( &
4046	( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) &
4047	* rho_air_zw(k) &
4048	- ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) &
4049	* rho_air_zw(k-1) &
4050	) * ddzw(k) * drho_air(k) &
4051	) * flag * dsig_diss &
4052	- c_2 * diss(k,j,i)**2 &
4053	/ ( e(k,j,i) + 1.0E-20_wp ) * flag
4054
4055	ENDDO
4056	ENDDO
4057	ENDDO
4058
4059	END SUBROUTINE diffusion_diss
4060
4061
4062	!------------------------------------------------------------------------------!
4063	! Description:
4064	! ------------
4065	!> Diffusion term for the TKE dissipation rate
4066	!> Cache-optimized version
4067	!------------------------------------------------------------------------------!
4068	SUBROUTINE diffusion_diss_ij( i, j )
4069
4070	USE arrays_3d, &
4071	ONLY: ddzu, ddzw, drho_air, rho_air_zw
4072
4073	USE grid_variables, &
4074	ONLY: ddx2, ddy2
4075
4076	IMPLICIT NONE
4077
4078	INTEGER(iwp) :: i !< running index x direction
4079	INTEGER(iwp) :: j !< running index y direction
4080	INTEGER(iwp) :: k !< running index z direction
4081
4082	REAL(wp) :: flag !< flag to mask topography
4083
4084	!
4085	!-- Calculate the mixing length (for dissipation)
4086	DO k = nzb+1, nzt
4087
4088	!
4089	!-- Predetermine flag to mask topography
4090	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4091
4092	!
4093	!-- Calculate the tendency term
4094	tend(k,j,i) = tend(k,j,i) + &
4095	( ( &
4096	( km(k,j,i)+km(k,j,i+1) ) * ( diss(k,j,i+1)-diss(k,j,i) ) &
4097	- ( km(k,j,i)+km(k,j,i-1) ) * ( diss(k,j,i)-diss(k,j,i-1) ) &
4098	) * ddx2 &
4099	+ ( &
4100	( km(k,j,i)+km(k,j+1,i) ) * ( diss(k,j+1,i)-diss(k,j,i) ) &
4101	- ( km(k,j,i)+km(k,j-1,i) ) * ( diss(k,j,i)-diss(k,j-1,i) ) &
4102	) * ddy2 &
4103	+ ( &
4104	( km(k,j,i)+km(k+1,j,i) ) * ( diss(k+1,j,i)-diss(k,j,i) ) * ddzu(k+1) &
4105	* rho_air_zw(k) &
4106	- ( km(k,j,i)+km(k-1,j,i) ) * ( diss(k,j,i)-diss(k-1,j,i) ) * ddzu(k) &
4107	* rho_air_zw(k-1) &
4108	) * ddzw(k) * drho_air(k) &
4109	) * flag * dsig_diss &
4110	- c_2 * diss(k,j,i)*2 / ( e(k,j,i) + 1.0E-20_wp ) flag
4111
4112	ENDDO
4113
4114	END SUBROUTINE diffusion_diss_ij
4115
4116
4117	!------------------------------------------------------------------------------!
4118	! Description:
4119	! ------------
4120	!> Calculate mixing length for LES mode.
4121	!------------------------------------------------------------------------------!
4122	SUBROUTINE mixing_length_les( i, j, k, l, ll, var, var_reference )
4123
4124	USE arrays_3d, &
4125	ONLY: dd2zu
4126
4127	USE control_parameters, &
4128	ONLY: atmos_ocean_sign, use_single_reference_value, &
4129	wall_adjustment, wall_adjustment_factor
4130
4131	IMPLICIT NONE
4132
4133	INTEGER(iwp) :: i !< loop index
4134	INTEGER(iwp) :: j !< loop index
4135	INTEGER(iwp) :: k !< loop index
4136
4137	REAL(wp) :: dvar_dz !< vertical gradient of var
4138	REAL(wp) :: l !< mixing length
4139	REAL(wp) :: l_stable !< mixing length according to stratification
4140	REAL(wp) :: ll !< adjusted l_grid
4141	REAL(wp) :: var_reference !< var at reference height
4142
4143	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4144
4145
4146	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4147	IF ( dvar_dz > 0.0_wp ) THEN
4148	IF ( use_single_reference_value ) THEN
4149	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
4150	/ SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
4151	ELSE
4152	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
4153	/ SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
4154	ENDIF
4155	ELSE
4156	l_stable = l_grid(k)
4157	ENDIF
4158	!
4159	!-- Adjustment of the mixing length
4160	IF ( wall_adjustment ) THEN
4161	l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
4162	ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
4163	ELSE
4164	l = MIN( l_grid(k), l_stable )
4165	ll = l_grid(k)
4166	ENDIF
4167
4168	END SUBROUTINE mixing_length_les
4169
4170
4171	!------------------------------------------------------------------------------!
4172	! Description:
4173	! ------------
4174	!> Calculate mixing length for RANS mode.
4175	!------------------------------------------------------------------------------!
4176	SUBROUTINE mixing_length_rans( i, j, k, l, l_diss, var, var_reference )
4177
4178	USE arrays_3d, &
4179	ONLY: dd2zu
4180
4181	USE control_parameters, &
4182	ONLY: atmos_ocean_sign, use_single_reference_value
4183
4184	IMPLICIT NONE
4185
4186	INTEGER(iwp) :: i !< loop index
4187	INTEGER(iwp) :: j !< loop index
4188	INTEGER(iwp) :: k !< loop index
4189
4190	REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector
4191	REAL(wp) :: dvar_dz !< vertical gradient of var
4192	REAL(wp) :: l !< mixing length
4193	REAL(wp) :: l_diss !< mixing length for dissipation
4194	REAL(wp) :: rif !< Richardson flux number
4195	REAL(wp) :: var_reference !< var at reference height
4196
4197	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4198
4199
4200	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4201
4202	duv2_dz2 = ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2 &
4203	+ ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2 &
4204	+ 1E-30_wp
4205
4206	IF ( use_single_reference_value ) THEN
4207	rif = g / var_reference * dvar_dz / duv2_dz2
4208	ELSE
4209	rif = g / var(k,j,i) * dvar_dz / duv2_dz2
4210	ENDIF
4211
4212	rif = MAX( rif, -5.0_wp )
4213	rif = MIN( rif, 1.0_wp )
4214
4215	!
4216	!-- Calculate diabatic mixing length using Dyer-profile functions
4217	IF ( rif >= 0.0_wp ) THEN
4218	l = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
4219	l_diss = l
4220	ELSE
4221	!
4222	!-- In case of unstable stratification, use mixing length of neutral case
4223	!-- for l, but consider profile functions for l_diss
4224	l = l_wall(k,j,i)
4225	l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
4226	ENDIF
4227
4228	END SUBROUTINE mixing_length_rans
4229
4230
4231	!------------------------------------------------------------------------------!
4232	! Description:
4233	! ------------
4234	!> Computation of the turbulent diffusion coefficients for momentum and heat.
4235	!> @bug unstable stratification is not properly considered for kh in rans mode.
4236	!------------------------------------------------------------------------------!
4237	SUBROUTINE tcm_diffusivities( var, var_reference )
4238
4239	USE control_parameters, &
4240	ONLY: bc_radiation_l, bc_radiation_n, bc_radiation_r, bc_radiation_s, &
4241	e_min
4242
4243	USE surface_layer_fluxes_mod, &
4244	ONLY: phi_m
4245
4246	USE surface_mod, &
4247	ONLY : bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v, &
4248	surf_usm_h, surf_usm_v
4249
4250	INTEGER(iwp) :: i !< loop index
4251	INTEGER(iwp) :: j !< loop index
4252	INTEGER(iwp) :: k !< loop index
4253	INTEGER(iwp) :: m !< loop index
4254	INTEGER(iwp) :: n !< loop index
4255
4256	REAL(wp) :: var_reference !< reference temperature
4257
4258	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4259
4260
4261	!
4262	!-- Introduce an optional minimum tke
4263	IF ( e_min > 0.0_wp ) THEN
4264	!$OMP PARALLEL DO PRIVATE(i,j,k)
4265	DO i = nxlg, nxrg
4266	DO j = nysg, nyng
4267	DO k = nzb+1, nzt
4268	e(k,j,i) = MAX( e(k,j,i), e_min ) * &
4269	MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4270	ENDDO
4271	ENDDO
4272	ENDDO
4273	ENDIF
4274
4275	!
4276	!-- Call default diffusivities routine. This is always used to calculate kh.
4277	CALL tcm_diffusivities_default( var, var_reference )
4278	!
4279	!-- Call dynamic subgrid model to calculate km.
4280	IF ( les_dynamic ) THEN
4281	CALL tcm_diffusivities_dynamic
4282	ENDIF
4283
4284	!
4285	!-- In RANS mode, use MOST to calculate km and kh within the surface layer.
4286	IF ( rans_tke_e ) THEN
4287	!
4288	!-- Upward facing surfaces
4289	!-- Default surfaces
4290	n = 0
4291	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4292	DO m = 1, surf_def_h(0)%ns
4293	i = surf_def_h(0)%i(m)
4294	j = surf_def_h(0)%j(m)
4295	k = surf_def_h(0)%k(m)
4296	km(k,j,i) = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) / &
4297	phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
4298	kh(k,j,i) = 1.35_wp * km(k,j,i)
4299	ENDDO
4300	!
4301	!-- Natural surfaces
4302	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4303	DO m = 1, surf_lsm_h%ns
4304	i = surf_lsm_h%i(m)
4305	j = surf_lsm_h%j(m)
4306	k = surf_lsm_h%k(m)
4307	km(k,j,i) = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) / &
4308	phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
4309	kh(k,j,i) = 1.35_wp * km(k,j,i)
4310	ENDDO
4311	!
4312	!-- Urban surfaces
4313	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4314	DO m = 1, surf_usm_h%ns
4315	i = surf_usm_h%i(m)
4316	j = surf_usm_h%j(m)
4317	k = surf_usm_h%k(m)
4318	km(k,j,i) = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) / &
4319	phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
4320	kh(k,j,i) = 1.35_wp * km(k,j,i)
4321	ENDDO
4322
4323	!
4324	!-- North-, south-, west and eastward facing surfaces
4325	!-- Do not consider stratification at these surfaces.
4326	DO n = 0, 3
4327	!
4328	!-- Default surfaces
4329	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4330	DO m = 1, surf_def_v(n)%ns
4331	i = surf_def_v(n)%i(m)
4332	j = surf_def_v(n)%j(m)
4333	k = surf_def_v(n)%k(m)
4334	km(k,j,i) = kappa * surf_def_v(n)%us(m) * surf_def_v(n)%z_mo(m)
4335	kh(k,j,i) = 1.35_wp * km(k,j,i)
4336	ENDDO
4337	!
4338	!-- Natural surfaces
4339	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4340	DO m = 1, surf_lsm_v(n)%ns
4341	i = surf_lsm_v(n)%i(m)
4342	j = surf_lsm_v(n)%j(m)
4343	k = surf_lsm_v(n)%k(m)
4344	km(k,j,i) = kappa * surf_lsm_v(n)%us(m) * surf_lsm_v(n)%z_mo(m)
4345	kh(k,j,i) = 1.35_wp * km(k,j,i)
4346	ENDDO
4347	!
4348	!-- Urban surfaces
4349	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
4350	DO m = 1, surf_usm_v(n)%ns
4351	i = surf_usm_v(n)%i(m)
4352	j = surf_usm_v(n)%j(m)
4353	k = surf_usm_v(n)%k(m)
4354	km(k,j,i) = kappa * surf_usm_v(n)%us(m) * surf_usm_v(n)%z_mo(m)
4355	kh(k,j,i) = 1.35_wp * km(k,j,i)
4356	ENDDO
4357	ENDDO
4358
4359	CALL exchange_horiz( km, nbgp )
4360	CALL exchange_horiz( kh, nbgp )
4361
4362	ENDIF
4363	!
4364	!-- Set boundary values (Neumann conditions)
4365	!-- Downward facing surfaces
4366	!$OMP PARALLEL DO PRIVATE(i,j,k)
4367	!$ACC PARALLEL LOOP PRIVATE(i,j,k) &
4368	!$ACC PRESENT(bc_h(1), kh, km)
4369	DO m = 1, bc_h(1)%ns
4370	i = bc_h(1)%i(m)
4371	j = bc_h(1)%j(m)
4372	k = bc_h(1)%k(m)
4373	km(k+1,j,i) = km(k,j,i)
4374	kh(k+1,j,i) = kh(k,j,i)
4375	ENDDO
4376	!
4377	!-- Downward facing surfaces
4378	!$OMP PARALLEL DO PRIVATE(i,j,k)
4379	!$ACC PARALLEL LOOP PRIVATE(i,j,k) &
4380	!$ACC PRESENT(bc_h(0), kh, km)
4381	DO m = 1, bc_h(0)%ns
4382	i = bc_h(0)%i(m)
4383	j = bc_h(0)%j(m)
4384	k = bc_h(0)%k(m)
4385	km(k-1,j,i) = km(k,j,i)
4386	kh(k-1,j,i) = kh(k,j,i)
4387	ENDDO
4388	!
4389	!-- Model top
4390	!$OMP PARALLEL DO
4391	!$ACC PARALLEL LOOP COLLAPSE(2) &
4392	!$ACC PRESENT(kh, km)
4393	DO i = nxlg, nxrg
4394	DO j = nysg, nyng
4395	km(nzt+1,j,i) = km(nzt,j,i)
4396	kh(nzt+1,j,i) = kh(nzt,j,i)
4397	ENDDO
4398	ENDDO
4399
4400	!
4401	!-- Set Neumann boundary conditions at the outflow boundaries in case of
4402	!-- non-cyclic lateral boundaries
4403	IF ( bc_radiation_l ) THEN
4404	km(:,:,nxl-1) = km(:,:,nxl)
4405	kh(:,:,nxl-1) = kh(:,:,nxl)
4406	ENDIF
4407	IF ( bc_radiation_r ) THEN
4408	km(:,:,nxr+1) = km(:,:,nxr)
4409	kh(:,:,nxr+1) = kh(:,:,nxr)
4410	ENDIF
4411	IF ( bc_radiation_s ) THEN
4412	km(:,nys-1,:) = km(:,nys,:)
4413	kh(:,nys-1,:) = kh(:,nys,:)
4414	ENDIF
4415	IF ( bc_radiation_n ) THEN
4416	km(:,nyn+1,:) = km(:,nyn,:)
4417	kh(:,nyn+1,:) = kh(:,nyn,:)
4418	ENDIF
4419
4420	END SUBROUTINE tcm_diffusivities
4421
4422
4423	!------------------------------------------------------------------------------!
4424	! Description:
4425	! ------------
4426	!> Computation of the turbulent diffusion coefficients for momentum and heat
4427	!> according to Prandtl-Kolmogorov.
4428	!------------------------------------------------------------------------------!
4429	SUBROUTINE tcm_diffusivities_default( var, var_reference )
4430
4431	USE arrays_3d, &
4432	ONLY: dd2zu
4433
4434	USE control_parameters, &
4435	ONLY: atmos_ocean_sign, use_single_reference_value, &
4436	wall_adjustment, wall_adjustment_factor
4437
4438	USE statistics, &
4439	ONLY : rmask, sums_l_l
4440
4441	IMPLICIT NONE
4442
4443	INTEGER(iwp) :: i !< loop index
4444	INTEGER(iwp) :: j !< loop index
4445	INTEGER(iwp) :: k !< loop index
4446	!$ INTEGER(iwp) :: omp_get_thread_num !< opemmp function to get thread number
4447	INTEGER(iwp) :: sr !< statistic region
4448	INTEGER(iwp) :: tn !< thread number
4449
4450	REAL(wp) :: flag !< topography flag
4451	REAL(wp) :: l !< mixing length
4452	REAL(wp) :: ll !< adjusted mixing length
4453	REAL(wp) :: var_reference !< reference temperature
4454	#ifdef _OPENACC
4455	!
4456	!-- From mixing_length_les:
4457	REAL(wp) :: l_stable !< mixing length according to stratification
4458	REAL(wp) :: dvar_dz !< vertical gradient of var
4459	#endif
4460
4461	REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature
4462
4463
4464	!
4465	!-- Default thread number in case of one thread
4466	tn = 0
4467
4468	!
4469	!-- Initialization for calculation of the mixing length profile
4470	!$ACC KERNELS PRESENT(sums_l_l)
4471	sums_l_l = 0.0_wp
4472	!$ACC END KERNELS
4473
4474	!
4475	!-- Compute the turbulent diffusion coefficient for momentum
4476	!$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
4477	!$ tn = omp_get_thread_num()
4478
4479	IF ( les_dynamic .OR. les_mw ) THEN
4480	!$OMP DO
4481	!$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(sr) &
4482	!$ACC PRIVATE(flag, dvar_dz, l_stable, l, ll) &
4483	!$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) &
4484	!$ACC PRESENT(kh, km, sums_l_l)
4485	DO i = nxlg, nxrg
4486	DO j = nysg, nyng
4487	DO k = nzb+1, nzt
4488
4489	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4490
4491	!
4492	!-- Determine the mixing length for LES closure
4493	#ifdef _OPENACC
4494	!
4495	!-- Cannot call subroutine mixing_length_les because after adding all required
4496	!-- OpenACC directives the execution crashed reliably when inside the called
4497	!-- subroutine. I'm not sure why that is...
4498	dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
4499	IF ( dvar_dz > 0.0_wp ) THEN
4500	IF ( use_single_reference_value ) THEN
4501	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
4502	/ SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
4503	ELSE
4504	l_stable = 0.76_wp * SQRT( e(k,j,i) ) &
4505	/ SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
4506	ENDIF
4507	ELSE
4508	l_stable = l_grid(k)
4509	ENDIF
4510	!
4511	!-- Adjustment of the mixing length
4512	IF ( wall_adjustment ) THEN
4513	l = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
4514	ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
4515	ELSE
4516	l = MIN( l_grid(k), l_stable )
4517	ll = l_grid(k)
4518	ENDIF
4519	#else
4520	CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
4521	#endif
4522
4523	!
4524	!-- Compute diffusion coefficients for momentum and heat
4525	km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
4526	kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
4527	!
4528	!-- Summation for averaged profile (cf. flow_statistics)
4529	DO sr = 0, statistic_regions
4530	sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) &
4531	* flag
4532	ENDDO
4533
4534	ENDDO
4535	ENDDO
4536	ENDDO
4537
4538	ELSEIF ( rans_tke_l ) THEN
4539
4540	!$OMP DO
4541	DO i = nxlg, nxrg
4542	DO j = nysg, nyng
4543	DO k = nzb+1, nzt
4544
4545	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4546	!
4547	!-- Mixing length for RANS mode with TKE-l closure
4548	CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
4549	!
4550	!-- Compute diffusion coefficients for momentum and heat
4551	km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
4552	kh(k,j,i) = km(k,j,i) / prandtl_number * flag
4553	!
4554	!-- Summation for averaged profile (cf. flow_statistics)
4555	DO sr = 0, statistic_regions
4556	sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr) &
4557	* flag
4558	ENDDO
4559
4560	ENDDO
4561	ENDDO
4562	ENDDO
4563
4564	ELSEIF ( rans_tke_e ) THEN
4565
4566	!$OMP DO
4567	DO i = nxlg, nxrg
4568	DO j = nysg, nyng
4569	DO k = nzb+1, nzt
4570
4571	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4572	!
4573	!-- Compute diffusion coefficients for momentum and heat
4574	km(k,j,i) = c_0*4 e(k,j,i)*2 / ( diss(k,j,i) + 1.0E-30_wp ) flag
4575	kh(k,j,i) = km(k,j,i) / prandtl_number * flag
4576	!
4577	!-- Summation for averaged profile of mixing length (cf. flow_statistics)
4578	DO sr = 0, statistic_regions
4579	sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + &
4580	c_0*3 e(k,j,i) * SQRT(e(k,j,i)) / &
4581	( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
4582	ENDDO
4583
4584	ENDDO
4585	ENDDO
4586	ENDDO
4587
4588	ENDIF
4589
4590	!$ACC KERNELS PRESENT(sums_l_l)
4591	sums_l_l(nzt+1,:,tn) = sums_l_l(nzt,:,tn) ! quasi boundary-condition for
4592	! data output
4593	!$ACC END KERNELS
4594	!$OMP END PARALLEL
4595
4596	END SUBROUTINE tcm_diffusivities_default
4597
4598
4599	!------------------------------------------------------------------------------!
4600	! Description:
4601	! ------------
4602	!> Calculates the eddy viscosity dynamically using the linear dynamic model
4603	!> according to
4604	!> Heinz, Stefan. "Realizability of dynamic subgrid-scale stress models via
4605	!> stochastic analysis."
4606	!> Monte Carlo Methods and Applications 14.4 (2008): 311-329.
4607	!>
4608	!> Furthermore dynamic bounds are used to limit the absolute value of c* as
4609	!> described in
4610	!> Mokhtarpoor, Reza, and Stefan Heinz. "Dynamic large eddy simulation:
4611	!> Stability via realizability." Physics of Fluids 29.10 (2017): 105104.
4612	!>
4613	!> @author Hauke Wurps
4614	!> @author BjÃ¶rn Maronga
4615	!------------------------------------------------------------------------------!
4616	SUBROUTINE tcm_diffusivities_dynamic
4617
4618	USE arrays_3d, &
4619	ONLY: ddzw, dzw, dd2zu, w, ug, vg
4620
4621	USE grid_variables, &
4622	ONLY : ddx, ddy, dx, dy
4623
4624	IMPLICIT NONE
4625
4626	INTEGER(iwp) :: i !< running index x-direction
4627	INTEGER(iwp) :: j !< running index y-direction
4628	INTEGER(iwp) :: k !< running index z-direction
4629	INTEGER(iwp) :: l !< running index
4630	INTEGER(iwp) :: m !< running index
4631
4632	REAL(wp) :: dudx !< Gradient of u-component in x-direction
4633	REAL(wp) :: dudy !< Gradient of u-component in y-direction
4634	REAL(wp) :: dudz !< Gradient of u-component in z-direction
4635	REAL(wp) :: dvdx !< Gradient of v-component in x-direction
4636	REAL(wp) :: dvdy !< Gradient of v-component in y-direction
4637	REAL(wp) :: dvdz !< Gradient of v-component in z-direction
4638	REAL(wp) :: dwdx !< Gradient of w-component in x-direction
4639	REAL(wp) :: dwdy !< Gradient of w-component in y-direction
4640	REAL(wp) :: dwdz !< Gradient of w-component in z-direction
4641
4642	REAL(wp) :: flag !< topography flag
4643
4644	REAL(wp) :: uc(-1:1,-1:1) !< u on grid center
4645	REAL(wp) :: vc(-1:1,-1:1) !< v on grid center
4646	REAL(wp) :: wc(-1:1,-1:1) !< w on grid center
4647
4648	REAL(wp) :: ut(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< test filtered u
4649	REAL(wp) :: vt(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< test filtered v
4650	REAL(wp) :: wt(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< test filtered w
4651
4652	REAL(wp) :: uct !< test filtered u on grid center
4653	REAL(wp) :: vct !< test filtered v on grid center
4654	REAL(wp) :: wct !< test filtered w on grid center
4655	REAL(wp) :: u2t !< test filtered u**2 on grid center
4656	REAL(wp) :: v2t !< test filtered v**2 on grid center
4657	REAL(wp) :: w2t !< test filtered w**2 on grid center
4658	REAL(wp) :: uvt !< test filtered u*v on grid center
4659	REAL(wp) :: uwt !< test filtered u*w on grid center
4660	REAL(wp) :: vwt !< test filtered v*w on grid center
4661
4662	REAL(wp) :: sd11 !< deviatoric shear tensor
4663	REAL(wp) :: sd22 !< deviatoric shear tensor
4664	REAL(wp) :: sd33 !<f deviatoric shear tensor
4665	REAL(wp) :: sd12 !< deviatoric shear tensor
4666	REAL(wp) :: sd13 !< deviatoric shear tensor
4667	REAL(wp) :: sd23 !< deviatoric shear tensor
4668
4669	REAL(wp) :: sd2 !< sum: sd_ij*sd_ij
4670
4671	REAL(wp) :: sdt11 !< filtered deviatoric shear tensor
4672	REAL(wp) :: sdt22 !< filtered deviatoric shear tensor
4673	REAL(wp) :: sdt33 !< filtered deviatoric shear tensor
4674	REAL(wp) :: sdt12 !< filtered deviatoric shear tensor
4675	REAL(wp) :: sdt13 !< filtered deviatoric shear tensor
4676	REAL(wp) :: sdt23 !< filtered deviatoric shear tensor
4677
4678	REAL(wp) :: sdt2 !< sum: sdt_ij*sdt_ij
4679
4680	REAL(wp) :: ld11 !< deviatoric stress tensor
4681	REAL(wp) :: ld22 !< deviatoric stress tensor
4682	REAL(wp) :: ld33 !< deviatoric stress tensor
4683	REAL(wp) :: ld12 !< deviatoric stress tensor
4684	REAL(wp) :: ld13 !< deviatoric stress tensor
4685	REAL(wp) :: ld23 !< deviatoric stress tensor
4686
4687	REAL(wp) :: lnn !< sum ld_nn
4688	REAL(wp) :: ldsd !< sum: ld_ij*sd_ij
4689
4690	REAL(wp) :: delta !< grid size
4691	REAL(wp) :: cst !< c*
4692	REAL(wp) :: cstnust_t !< product cnu
4693	REAL(wp) :: cst_max !< bounds of c*
4694
4695	REAL(wp), PARAMETER :: fac_cmax = 23.0_wp/(24.0_wp*sqrt(3.0_wp)) !< constant
4696
4697	!
4698	!-- velocities on grid centers:
4699	CALL tcm_box_filter_2d( u, ut )
4700	CALL tcm_box_filter_2d( v, vt )
4701	CALL tcm_box_filter_2d( w, wt )
4702
4703	DO i = nxl, nxr
4704	DO j = nys, nyn
4705	DO k = nzb+1, nzt
4706
4707	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
4708
4709	!
4710	!-- Compute the deviatoric shear tensor s_ij on grid centers:
4711	!-- s_ij = 0.5 * ( du_i/dx_j + du_j/dx_i )
4712	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
4713	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
4714	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
4715	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
4716	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
4717
4718	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
4719	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
4720	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
4721	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
4722	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
4723
4724	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
4725	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
4726	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
4727	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
4728	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
4729
4730	sd11 = dudx
4731	sd22 = dvdy
4732	sd33 = dwdz
4733	sd12 = 0.5_wp * ( dudy + dvdx )
4734	sd13 = 0.5_wp * ( dudz + dwdx )
4735	sd23 = 0.5_wp * ( dvdz + dwdy )
4736	!
4737	!-- sum: sd_ij*sd_ij
4738	sd2 = sd112 + sd222 + sd33**2 &
4739	+ 2.0_wp * ( sd122 + sd132 + sd23**2 )
4740	!
4741	!-- The filtered velocities are needed to calculate the filtered shear
4742	!-- tensor: sdt_ij = 0.5 * ( dut_i/dx_j + dut_j/dx_i )
4743	dudx = ( ut(k,j,i+1) - ut(k,j,i) ) * ddx
4744	dudy = 0.25_wp * ( ut(k,j+1,i) + ut(k,j+1,i+1) - &
4745	ut(k,j-1,i) - ut(k,j-1,i+1) ) * ddy
4746	dudz = 0.5_wp * ( ut(k+1,j,i) + ut(k+1,j,i+1) - &
4747	ut(k-1,j,i) - ut(k-1,j,i+1) ) * dd2zu(k)
4748
4749	dvdx = 0.25_wp * ( vt(k,j,i+1) + vt(k,j+1,i+1) - &
4750	vt(k,j,i-1) - vt(k,j+1,i-1) ) * ddx
4751	dvdy = ( vt(k,j+1,i) - vt(k,j,i) ) * ddy
4752	dvdz = 0.5_wp * ( vt(k+1,j,i) + vt(k+1,j+1,i) - &
4753	vt(k-1,j,i) - vt(k-1,j+1,i) ) * dd2zu(k)
4754
4755	dwdx = 0.25_wp * ( wt(k,j,i+1) + wt(k-1,j,i+1) - &
4756	wt(k,j,i-1) - wt(k-1,j,i-1) ) * ddx
4757	dwdy = 0.25_wp * ( wt(k,j+1,i) + wt(k-1,j+1,i) - &
4758	wt(k,j-1,i) - wt(k-1,j-1,i) ) * ddy
4759	dwdz = ( wt(k,j,i) - wt(k-1,j,i) ) * ddzw(k)
4760
4761	sdt11 = dudx
4762	sdt22 = dvdy
4763	sdt33 = dwdz
4764	sdt12 = 0.5_wp * ( dudy + dvdx )
4765	sdt13 = 0.5_wp * ( dudz + dwdx )
4766	sdt23 = 0.5_wp * ( dvdz + dwdy )
4767	!
4768	!-- sum: sd_ij*sd_ij
4769	sdt2 = sdt112 + sdt222 + sdt33**2 &
4770	+ 2.0_wp * ( sdt122 + sdt132 + sdt23**2 )
4771	!
4772	!-- Need filtered velocities and filtered squared velocities on grid
4773	!-- centers. Substraction of geostrophic velocity helps to avoid
4774	!-- numerical errors in the expression <u*2> - <u><u>, which can be
4775	!-- very small (<...> means filtered).
4776	DO l = -1, 1
4777	DO m = -1, 1
4778	uc(l,m) = 0.5_wp * ( u(k,j+l,i+m) + u(k,j+l,i+m+1) ) - ug(k)
4779	vc(l,m) = 0.5_wp * ( v(k,j+l,i+m) + v(k,j+l+1,i+m) ) - vg(k)
4780	wc(l,m) = 0.5_wp * ( w(k-1,j+l,i+m) + w(k,j+l,i+m) )
4781	ENDDO
4782	ENDDO
4783
4784	CALL tcm_box_filter_2d( uc, uct )
4785	CALL tcm_box_filter_2d( vc, vct )
4786	CALL tcm_box_filter_2d( wc, wct )
4787	CALL tcm_box_filter_2d( uc**2, u2t )
4788	CALL tcm_box_filter_2d( vc**2, v2t )
4789	CALL tcm_box_filter_2d( wc**2, w2t )
4790	CALL tcm_box_filter_2d( uc*vc, uvt )
4791	CALL tcm_box_filter_2d( uc*wc, uwt )
4792	CALL tcm_box_filter_2d( vc*wc, vwt )
4793
4794	ld11 = u2t - uct*uct
4795	ld22 = v2t - vct*vct
4796	ld33 = w2t - wct*wct
4797	ld12 = uvt - uct*vct
4798	ld13 = uwt - uct*wct
4799	ld23 = vwt - vct*wct
4800
4801	lnn = ld11 + ld22 + ld33
4802	!
4803	!-- Substract trace to get deviatoric resolved stress
4804	ld11 = ld11 - lnn / 3.0_wp
4805	ld22 = ld22 - lnn / 3.0_wp
4806	ld33 = ld33 - lnn / 3.0_wp
4807
4808	ldsd = ld11 * sdt11 + ld22 * sdt22 + ld33 * sdt33 + &
4809	2.0_wp * ( ld12 * sdt12 + ld13 * sdt13 + ld23 * sdt23 )
4810	!
4811	!-- c* nu*^T is SGS viscosity on test filter level:
4812	cstnust_t = -ldsd / ( sdt2 + 1.0E-20_wp )
4813	!
4814	!-- The model was only tested for an isotropic grid. The following
4815	!-- expression was a recommendation of Stefan Heinz.
4816	delta = MAX( dx, dy, dzw(k) )
4817
4818	IF ( lnn <= 0.0_wp ) THEN
4819	cst = 0.0_wp
4820	ELSE
4821	cst = cstnust_t / &
4822	( 4.0_wp * delta * SQRT( lnn / 2.0_wp ) + 1.0E-20_wp )
4823	ENDIF
4824
4825	!
4826	!-- Calculate border according to Mokhtarpoor and Heinz (2017)
4827	cst_max = fac_cmax * SQRT( e(k,j,i) ) / &
4828	( delta * SQRT( 2.0_wp * sd2 ) + 1.0E-20_wp )
4829
4830	IF ( ABS( cst ) > cst_max ) THEN
4831	cst = cst_max * cst / ABS( cst )
4832	ENDIF
4833
4834	km(k,j,i) = cst * delta * SQRT( e(k,j,i) ) * flag
4835
4836	ENDDO
4837	ENDDO
4838	ENDDO
4839
4840	END SUBROUTINE tcm_diffusivities_dynamic
4841
4842
4843	!------------------------------------------------------------------------------!
4844	! Description:
4845	! ------------
4846	!> This subroutine acts as a box filter with filter width 2 * dx.
4847	!> Output is only one point.
4848	!------------------------------------------------------------------------------!
4849	SUBROUTINE tcm_box_filter_2d_single( var, var_fil )
4850
4851	IMPLICIT NONE
4852
4853	REAL(wp) :: var(-1:1,-1:1) !< variable to be filtered
4854	REAL(wp) :: var_fil !< filtered variable
4855	!
4856	!-- It is assumed that a box with a side length of 2 * dx and centered at the
4857	!-- variable*s position contains one half of the four closest neigbours and one
4858	!-- forth of the diagonally closest neighbours.
4859	var_fil = 0.25_wp * ( var(0,0) + &
4860	0.5_wp * ( var(0,1) + var(1,0) + &
4861	var(0,-1) + var(-1,0) ) + &
4862	0.25_wp * ( var(1,1) + var(1,-1) + &
4863	var(-1,1) + var(-1,-1) ) )
4864
4865	END SUBROUTINE tcm_box_filter_2d_single
4866
4867	!------------------------------------------------------------------------------!
4868	! Description:
4869	! ------------
4870	!> This subroutine acts as a box filter with filter width 2 * dx.
4871	!> The filtered variable var_fil is on the same grid as var.
4872	!------------------------------------------------------------------------------!
4873	SUBROUTINE tcm_box_filter_2d_array( var, var_fil )
4874
4875	IMPLICIT NONE
4876
4877	INTEGER(iwp) :: i !< running index x-direction
4878	INTEGER(iwp) :: j !< running index y-direction
4879	INTEGER(iwp) :: k !< running index z-direction
4880
4881	REAL(wp) :: var(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< variable to be filtered
4882	REAL(wp) :: var_fil(nzb:nzt+1,nysg:nyng,nxlg:nxrg) !< filtered variable
4883	!
4884	!-- It is assumed that a box with a side length of 2 * dx and centered at the
4885	!-- variable's position contains one half of the four closest neigbours and one
4886	!-- forth of the diagonally closest neighbours.
4887	DO i = nxlg+1, nxrg-1
4888	DO j = nysg+1, nyng-1
4889	DO k = nzb, nzt+1
4890	var_fil(k,j,i) = 0.25_wp * ( var(k,j,i) + &
4891	0.5_wp * ( var(k,j,i+1) + var(k,j+1,i) + &
4892	var(k,j,i-1) + var(k,j-1,i) ) +&
4893	0.25_wp * ( var(k,j+1,i+1) + var(k,j+1,i-1) + &
4894	var(k,j-1,i+1) + var(k,j-1,i-1) ) )
4895	END DO
4896	END DO
4897	END DO
4898
4899	END SUBROUTINE tcm_box_filter_2d_array
4900
4901
4902	!------------------------------------------------------------------------------!
4903	! Description:
4904	! ------------
4905	!> Swapping of timelevels.
4906	!------------------------------------------------------------------------------!
4907	SUBROUTINE tcm_swap_timelevel ( mod_count )
4908
4909	IMPLICIT NONE
4910
4911
4912	INTEGER, INTENT(IN) :: mod_count !< flag defining where pointers point to
4913
4914
4915	SELECT CASE ( mod_count )
4916
4917	CASE ( 0 )
4918
4919	IF ( .NOT. constant_diffusion ) THEN
4920	e => e_1; e_p => e_2
4921	ENDIF
4922
4923	IF ( rans_tke_e ) THEN
4924	diss => diss_1; diss_p => diss_2
4925	ENDIF
4926
4927	CASE ( 1 )
4928
4929	IF ( .NOT. constant_diffusion ) THEN
4930	e => e_2; e_p => e_1
4931	ENDIF
4932
4933	IF ( rans_tke_e ) THEN
4934	diss => diss_2; diss_p => diss_1
4935	ENDIF
4936
4937	END SELECT
4938
4939	END SUBROUTINE tcm_swap_timelevel
4940
4941
4942	END MODULE turbulence_closure_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |