Home

Context Navigation

turbulence_closure_mod.f90 @ 3045

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