Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |