Home

Context Navigation

turbulence_closure_mod.f90 @ 3719

Last change on this file since 3719 was 3719, checked in by kanani, 6 years ago
Correct and clean-up cpu_logs, some overlapping counts (chemistry_model_mod, disturb_heatflux, large_scale_forcing_nudging_mod, ocean_mod, palm, prognostic_equations, synthetic_turbulence_generator_mod, time_integration, time_integration_spinup, turbulence_closure_mod)
Property svn:keywords set to `Id`
File size: 193.4 KB

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