Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |