Home

Context Navigation

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

Last change on this file since 3634 was 3634, checked in by knoop, 6 years ago
OpenACC port for SPEC
Property svn:keywords set to `Id`
File size: 195.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |