Home

Context Navigation

turbulence_closure_mod.f90 @ 3684

Last change on this file since 3684 was 3684, checked in by knoop, 6 years ago
Enabled module_interface_actions in time_integration and prognostic_equations
Property svn:keywords set to `Id`
File size: 193.3 KB

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