Home

Context Navigation

turbulence_closure_mod.f90 @ 4048

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