Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |