Home

Context Navigation

source: palm/trunk/SOURCE/production_e.f90 @ 1985

Last change on this file since 1985 was 1874, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 89.4 KB

Line
1	!> @file production_e.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2016 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	!
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: production_e.f90 1874 2016-04-18 14:51:10Z suehring $
26	!
27	! 1873 2016-04-18 14:50:06Z maronga
28	! Module renamed (removed _mod)
29	!
30	!
31	! 1850 2016-04-08 13:29:27Z maronga
32	! Module renamed
33	!
34	!
35	! 1691 2015-10-26 16:17:44Z maronga
36	! Renamed prandtl_layer to constant_flux_layer.
37	!
38	! 1682 2015-10-07 23:56:08Z knoop
39	! Code annotations made doxygen readable
40	!
41	! 1374 2014-04-25 12:55:07Z raasch
42	! nzb_s_outer removed from acc-present-list
43	!
44	! 1353 2014-04-08 15:21:23Z heinze
45	! REAL constants provided with KIND-attribute
46	!
47	! 1342 2014-03-26 17:04:47Z kanani
48	! REAL constants defined as wp-kind
49	!
50	! 1320 2014-03-20 08:40:49Z raasch
51	! ONLY-attribute added to USE-statements,
52	! kind-parameters added to all INTEGER and REAL declaration statements,
53	! kinds are defined in new module kinds,
54	! old module precision_kind is removed,
55	! revision history before 2012 removed,
56	! comment fields (!:) to be used for variable explanations added to
57	! all variable declaration statements
58	!
59	! 1257 2013-11-08 15:18:40Z raasch
60	! openacc loop and loop vector clauses removed, declare create moved after
61	! the FORTRAN declaration statement
62	!
63	! 1179 2013-06-14 05:57:58Z raasch
64	! use_reference renamed use_single_reference_value
65	!
66	! 1128 2013-04-12 06:19:32Z raasch
67	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
68	! j_north
69	!
70	! 1036 2012-10-22 13:43:42Z raasch
71	! code put under GPL (PALM 3.9)
72	!
73	! 1015 2012-09-27 09:23:24Z raasch
74	! accelerator version (*_acc) added
75	!
76	! 1007 2012-09-19 14:30:36Z franke
77	! Bugfix: calculation of buoyancy production has to consider the liquid water
78	! mixing ratio in case of cloud droplets
79	!
80	! 940 2012-07-09 14:31:00Z raasch
81	! TKE production by buoyancy can be switched off in case of runs with pure
82	! neutral stratification
83	!
84	! Revision 1.1 1997/09/19 07:45:35 raasch
85	! Initial revision
86	!
87	!
88	! Description:
89	! ------------
90	!> Production terms (shear + buoyancy) of the TKE.
91	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
92	!> not considered well!
93	!------------------------------------------------------------------------------!
94	MODULE production_e_mod
95
96
97	USE wall_fluxes_mod, &
98	ONLY: wall_fluxes_e, wall_fluxes_e_acc
99
100	USE kinds
101
102	PRIVATE
103	PUBLIC production_e, production_e_acc, production_e_init
104
105	LOGICAL, SAVE :: first_call = .TRUE. !<
106
107	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0 !<
108	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: v_0 !<
109
110	INTERFACE production_e
111	MODULE PROCEDURE production_e
112	MODULE PROCEDURE production_e_ij
113	END INTERFACE production_e
114
115	INTERFACE production_e_acc
116	MODULE PROCEDURE production_e_acc
117	END INTERFACE production_e_acc
118
119	INTERFACE production_e_init
120	MODULE PROCEDURE production_e_init
121	END INTERFACE production_e_init
122
123	CONTAINS
124
125
126	!------------------------------------------------------------------------------!
127	! Description:
128	! ------------
129	!> Call for all grid points
130	!------------------------------------------------------------------------------!
131	SUBROUTINE production_e
132
133	USE arrays_3d, &
134	ONLY: ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho, shf, &
135	tend, tswst, u, v, vpt, w
136
137	USE cloud_parameters, &
138	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
139
140	USE control_parameters, &
141	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
142	humidity, kappa, neutral, ocean, pt_reference, &
143	rho_reference, use_single_reference_value, &
144	use_surface_fluxes, use_top_fluxes
145
146	USE grid_variables, &
147	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
148
149	USE indices, &
150	ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, &
151	nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
152
153	IMPLICIT NONE
154
155	INTEGER(iwp) :: i !<
156	INTEGER(iwp) :: j !<
157	INTEGER(iwp) :: k !<
158
159	REAL(wp) :: def !<
160	REAL(wp) :: dudx !<
161	REAL(wp) :: dudy !<
162	REAL(wp) :: dudz !<
163	REAL(wp) :: dvdx !<
164	REAL(wp) :: dvdy !<
165	REAL(wp) :: dvdz !<
166	REAL(wp) :: dwdx !<
167	REAL(wp) :: dwdy !<
168	REAL(wp) :: dwdz !<
169	REAL(wp) :: k1 !<
170	REAL(wp) :: k2 !<
171	REAL(wp) :: km_neutral !<
172	REAL(wp) :: theta !<
173	REAL(wp) :: temp !<
174
175	! REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs
176	REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !<
177	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
178	REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !<
179	REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !<
180
181	!
182	!-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
183	!-- vertical walls, if neccessary
184	!-- So far, results are slightly different from the ij-Version.
185	!-- Therefore, ij-Version is called further below within the ij-loops.
186	! IF ( topography /= 'flat' ) THEN
187	! CALL wall_fluxes_e( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )
188	! CALL wall_fluxes_e( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )
189	! CALL wall_fluxes_e( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )
190	! CALL wall_fluxes_e( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )
191	! ENDIF
192
193
194	DO i = nxl, nxr
195
196	!
197	!-- Calculate TKE production by shear
198	DO j = nys, nyn
199	DO k = nzb_diff_s_outer(j,i), nzt
200
201	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
202	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
203	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
204	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
205	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
206
207	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
208	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
209	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
210	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
211	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
212
213	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
214	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
215	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
216	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
217	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
218
219	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
220	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
221	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
222
223	IF ( def < 0.0_wp ) def = 0.0_wp
224
225	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
226
227	ENDDO
228	ENDDO
229
230	IF ( constant_flux_layer ) THEN
231
232	!
233	!-- Position beneath wall
234	!-- (2) - Will allways be executed.
235	!-- 'bottom and wall: use u_0,v_0 and wall functions'
236	DO j = nys, nyn
237
238	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
239	THEN
240
241	k = nzb_diff_s_inner(j,i) - 1
242	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
243	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
244	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
245	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
246	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
247	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
248	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
249
250	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
251	!
252	!-- Inconsistency removed: as the thermal stratification is
253	!-- not taken into account for the evaluation of the wall
254	!-- fluxes at vertical walls, the eddy viscosity km must not
255	!-- be used for the evaluation of the velocity gradients dudy
256	!-- and dwdy
257	!-- Note: The validity of the new method has not yet been
258	!-- shown, as so far no suitable data for a validation
259	!-- has been available
260	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
261	usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
262	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
263	wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
264	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25_wp &
265	0.5_wp * dy
266	IF ( km_neutral > 0.0_wp ) THEN
267	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
268	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
269	ELSE
270	dudy = 0.0_wp
271	dwdy = 0.0_wp
272	ENDIF
273	ELSE
274	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
275	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
276	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
277	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
278	ENDIF
279
280	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
281	!
282	!-- Inconsistency removed: as the thermal stratification is
283	!-- not taken into account for the evaluation of the wall
284	!-- fluxes at vertical walls, the eddy viscosity km must not
285	!-- be used for the evaluation of the velocity gradients dvdx
286	!-- and dwdx
287	!-- Note: The validity of the new method has not yet been
288	!-- shown, as so far no suitable data for a validation
289	!-- has been available
290	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
291	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
292	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
293	wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
294	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25_wp &
295	0.5_wp * dx
296	IF ( km_neutral > 0.0_wp ) THEN
297	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
298	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
299	ELSE
300	dvdx = 0.0_wp
301	dwdx = 0.0_wp
302	ENDIF
303	ELSE
304	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
305	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
306	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
307	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
308	ENDIF
309
310	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
311	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
312	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
313
314	IF ( def < 0.0_wp ) def = 0.0_wp
315
316	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
317
318
319	!
320	!-- (3) - will be executed only, if there is at least one level
321	!-- between (2) and (4), i.e. the topography must have a
322	!-- minimum height of 2 dz. Wall fluxes for this case have
323	!-- already been calculated for (2).
324	!-- 'wall only: use wall functions'
325
326	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
327
328	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
329	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
330	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
331	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
332	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
333	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
334	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
335
336	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
337	!
338	!-- Inconsistency removed: as the thermal stratification
339	!-- is not taken into account for the evaluation of the
340	!-- wall fluxes at vertical walls, the eddy viscosity km
341	!-- must not be used for the evaluation of the velocity
342	!-- gradients dudy and dwdy
343	!-- Note: The validity of the new method has not yet
344	!-- been shown, as so far no suitable data for a
345	!-- validation has been available
346	km_neutral = kappa * ( usvs(k)**2 + &
347	wsvs(k)2 )0.25_wp * 0.5_wp * dy
348	IF ( km_neutral > 0.0_wp ) THEN
349	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
350	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
351	ELSE
352	dudy = 0.0_wp
353	dwdy = 0.0_wp
354	ENDIF
355	ELSE
356	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
357	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
358	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
359	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
360	ENDIF
361
362	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
363	!
364	!-- Inconsistency removed: as the thermal stratification
365	!-- is not taken into account for the evaluation of the
366	!-- wall fluxes at vertical walls, the eddy viscosity km
367	!-- must not be used for the evaluation of the velocity
368	!-- gradients dvdx and dwdx
369	!-- Note: The validity of the new method has not yet
370	!-- been shown, as so far no suitable data for a
371	!-- validation has been available
372	km_neutral = kappa * ( vsus(k)**2 + &
373	wsus(k)2 )0.25_wp * 0.5_wp * dx
374	IF ( km_neutral > 0.0_wp ) THEN
375	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
376	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
377	ELSE
378	dvdx = 0.0_wp
379	dwdx = 0.0_wp
380	ENDIF
381	ELSE
382	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
383	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
384	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
385	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
386	ENDIF
387
388	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
389	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
390	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
391
392	IF ( def < 0.0_wp ) def = 0.0_wp
393
394	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
395
396	ENDDO
397
398	ENDIF
399
400	ENDDO
401
402	!
403	!-- (4) - will allways be executed.
404	!-- 'special case: free atmosphere' (as for case (0))
405	DO j = nys, nyn
406
407	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
408	THEN
409
410	k = nzb_diff_s_outer(j,i)-1
411
412	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
413	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
414	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
415	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
416	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
417
418	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
419	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
420	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
421	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
422	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
423
424	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
425	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
426	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
427	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
428	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
429
430	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
431	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
432	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
433
434	IF ( def < 0.0_wp ) def = 0.0_wp
435
436	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
437
438	ENDIF
439
440	ENDDO
441
442	!
443	!-- Position without adjacent wall
444	!-- (1) - will allways be executed.
445	!-- 'bottom only: use u_0,v_0'
446	DO j = nys, nyn
447
448	IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
449	THEN
450
451	k = nzb_diff_s_inner(j,i)-1
452
453	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
454	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
455	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
456	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
457	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
458
459	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
460	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
461	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
462	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
463	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
464
465	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
466	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
467	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
468	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
469	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
470
471	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
472	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
473	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
474
475	IF ( def < 0.0_wp ) def = 0.0_wp
476
477	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
478
479	ENDIF
480
481	ENDDO
482
483	ELSEIF ( use_surface_fluxes ) THEN
484
485	DO j = nys, nyn
486
487	k = nzb_diff_s_outer(j,i)-1
488
489	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
490	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
491	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
492	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
493	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
494
495	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
496	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
497	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
498	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
499	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
500
501	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
502	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
503	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
504	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
505	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
506
507	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
508	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
509	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
510
511	IF ( def < 0.0_wp ) def = 0.0_wp
512
513	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
514
515	ENDDO
516
517	ENDIF
518
519	!
520	!-- If required, calculate TKE production by buoyancy
521	IF ( .NOT. neutral ) THEN
522
523	IF ( .NOT. humidity ) THEN
524
525	IF ( use_single_reference_value ) THEN
526
527	IF ( ocean ) THEN
528	!
529	!-- So far in the ocean no special treatment of density flux
530	!-- in the bottom and top surface layer
531	DO j = nys, nyn
532	DO k = nzb_s_inner(j,i)+1, nzt
533	tend(k,j,i) = tend(k,j,i) + &
534	kh(k,j,i) * g / rho_reference * &
535	( rho(k+1,j,i) - rho(k-1,j,i) ) * &
536	dd2zu(k)
537	ENDDO
538	ENDDO
539
540	ELSE
541
542	DO j = nys, nyn
543	DO k = nzb_diff_s_inner(j,i), nzt_diff
544	tend(k,j,i) = tend(k,j,i) - &
545	kh(k,j,i) * g / pt_reference * &
546	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
547	dd2zu(k)
548	ENDDO
549
550	IF ( use_surface_fluxes ) THEN
551	k = nzb_diff_s_inner(j,i)-1
552	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
553	shf(j,i)
554	ENDIF
555
556	IF ( use_top_fluxes ) THEN
557	k = nzt
558	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
559	tswst(j,i)
560	ENDIF
561	ENDDO
562
563	ENDIF
564
565	ELSE
566
567	IF ( ocean ) THEN
568	!
569	!-- So far in the ocean no special treatment of density flux
570	!-- in the bottom and top surface layer
571	DO j = nys, nyn
572	DO k = nzb_s_inner(j,i)+1, nzt
573	tend(k,j,i) = tend(k,j,i) + &
574	kh(k,j,i) * g / rho(k,j,i) * &
575	( rho(k+1,j,i) - rho(k-1,j,i) ) * &
576	dd2zu(k)
577	ENDDO
578	ENDDO
579
580	ELSE
581
582	DO j = nys, nyn
583	DO k = nzb_diff_s_inner(j,i), nzt_diff
584	tend(k,j,i) = tend(k,j,i) - &
585	kh(k,j,i) * g / pt(k,j,i) * &
586	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
587	dd2zu(k)
588	ENDDO
589
590	IF ( use_surface_fluxes ) THEN
591	k = nzb_diff_s_inner(j,i)-1
592	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
593	shf(j,i)
594	ENDIF
595
596	IF ( use_top_fluxes ) THEN
597	k = nzt
598	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
599	tswst(j,i)
600	ENDIF
601	ENDDO
602
603	ENDIF
604
605	ENDIF
606
607	ELSE
608
609	DO j = nys, nyn
610
611	DO k = nzb_diff_s_inner(j,i), nzt_diff
612
613	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
614	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
615	k2 = 0.61_wp * pt(k,j,i)
616	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
617	g / vpt(k,j,i) * &
618	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
619	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
620	) * dd2zu(k)
621	ELSE IF ( cloud_physics ) THEN
622	IF ( ql(k,j,i) == 0.0_wp ) THEN
623	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
624	k2 = 0.61_wp * pt(k,j,i)
625	ELSE
626	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
627	temp = theta * t_d_pt(k)
628	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
629	( q(k,j,i) - ql(k,j,i) ) * &
630	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
631	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
632	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
633	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
634	ENDIF
635	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
636	g / vpt(k,j,i) * &
637	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
638	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
639	) * dd2zu(k)
640	ELSE IF ( cloud_droplets ) THEN
641	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
642	k2 = 0.61_wp * pt(k,j,i)
643	tend(k,j,i) = tend(k,j,i) - &
644	kh(k,j,i) * g / vpt(k,j,i) * &
645	( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &
646	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
647	pt(k,j,i) * ( ql(k+1,j,i) - &
648	ql(k-1,j,i) ) ) * dd2zu(k)
649	ENDIF
650
651	ENDDO
652
653	ENDDO
654
655	IF ( use_surface_fluxes ) THEN
656
657	DO j = nys, nyn
658
659	k = nzb_diff_s_inner(j,i)-1
660
661	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
662	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
663	k2 = 0.61_wp * pt(k,j,i)
664	ELSE IF ( cloud_physics ) THEN
665	IF ( ql(k,j,i) == 0.0_wp ) THEN
666	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
667	k2 = 0.61_wp * pt(k,j,i)
668	ELSE
669	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
670	temp = theta * t_d_pt(k)
671	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
672	( q(k,j,i) - ql(k,j,i) ) * &
673	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
674	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
675	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
676	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
677	ENDIF
678	ELSE IF ( cloud_droplets ) THEN
679	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
680	k2 = 0.61_wp * pt(k,j,i)
681	ENDIF
682
683	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
684	( k1* shf(j,i) + k2 * qsws(j,i) )
685	ENDDO
686
687	ENDIF
688
689	IF ( use_top_fluxes ) THEN
690
691	DO j = nys, nyn
692
693	k = nzt
694
695	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
696	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
697	k2 = 0.61_wp * pt(k,j,i)
698	ELSE IF ( cloud_physics ) THEN
699	IF ( ql(k,j,i) == 0.0_wp ) THEN
700	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
701	k2 = 0.61_wp * pt(k,j,i)
702	ELSE
703	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
704	temp = theta * t_d_pt(k)
705	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
706	( q(k,j,i) - ql(k,j,i) ) * &
707	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
708	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
709	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
710	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
711	ENDIF
712	ELSE IF ( cloud_droplets ) THEN
713	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
714	k2 = 0.61_wp * pt(k,j,i)
715	ENDIF
716
717	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
718	( k1* tswst(j,i) + k2 * qswst(j,i) )
719	ENDDO
720
721	ENDIF
722
723	ENDIF
724
725	ENDIF
726
727	ENDDO
728
729	END SUBROUTINE production_e
730
731
732	!------------------------------------------------------------------------------!
733	! Description:
734	! ------------
735	!> Call for all grid points - accelerator version
736	!------------------------------------------------------------------------------!
737	SUBROUTINE production_e_acc
738
739	USE arrays_3d, &
740	ONLY: ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho, shf, &
741	tend, tswst, u, v, vpt, w
742
743	USE cloud_parameters, &
744	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
745
746	USE control_parameters, &
747	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
748	humidity, kappa, neutral, ocean, pt_reference, &
749	rho_reference, topography, use_single_reference_value, &
750	use_surface_fluxes, use_top_fluxes
751
752	USE grid_variables, &
753	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
754
755	USE indices, &
756	ONLY: i_left, i_right, j_north, j_south, nxl, nxr, nys, nyn, nzb, &
757	nzb_diff_s_inner, nzb_diff_s_outer, nzb_s_inner, nzt, &
758	nzt_diff
759
760	IMPLICIT NONE
761
762	INTEGER(iwp) :: i !<
763	INTEGER(iwp) :: j !<
764	INTEGER(iwp) :: k !<
765
766	REAL(wp) :: def !<
767	REAL(wp) :: dudx !<
768	REAL(wp) :: dudy !<
769	REAL(wp) :: dudz !<
770	REAL(wp) :: dvdx !<
771	REAL(wp) :: dvdy !<
772	REAL(wp) :: dvdz !<
773	REAL(wp) :: dwdx !<
774	REAL(wp) :: dwdy !<
775	REAL(wp) :: dwdz !<
776	REAL(wp) :: k1 !<
777	REAL(wp) :: k2 !<
778	REAL(wp) :: km_neutral !<
779	REAL(wp) :: theta !<
780	REAL(wp) :: temp !<
781
782	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs !<
783	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: vsus !<
784	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsus !<
785	REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: wsvs !<
786	!$acc declare create ( usvs, vsus, wsus, wsvs )
787
788	!
789	!-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
790	!-- vertical walls, if neccessary
791	!-- CAUTION: results are slightly different from the ij-version!!
792	!-- ij-version should be called further below within the ij-loops!!
793	IF ( topography /= 'flat' ) THEN
794	CALL wall_fluxes_e_acc( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )
795	CALL wall_fluxes_e_acc( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )
796	CALL wall_fluxes_e_acc( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )
797	CALL wall_fluxes_e_acc( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )
798	ENDIF
799
800
801	!
802	!-- Calculate TKE production by shear
803	!$acc kernels present( ddzw, dd2zu, kh, km, nzb_diff_s_inner, nzb_diff_s_outer ) &
804	!$acc present( nzb_s_inner, pt, q, ql, qsws, qswst, rho ) &
805	!$acc present( shf, tend, tswst, u, v, vpt, w, wall_e_x, wall_e_y ) &
806	!$acc copyin( u_0, v_0 )
807	DO i = i_left, i_right
808	DO j = j_south, j_north
809	DO k = 1, nzt
810
811	IF ( k >= nzb_diff_s_outer(j,i) ) THEN
812
813	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
814	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
815	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
816	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
817	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
818
819	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
820	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
821	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
822	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
823	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
824
825	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
826	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
827	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
828	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
829	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
830
831	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
832	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
833	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
834
835	IF ( def < 0.0_wp ) def = 0.0_wp
836
837	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
838
839	ENDIF
840
841	ENDDO
842	ENDDO
843	ENDDO
844
845	IF ( constant_flux_layer ) THEN
846
847	!
848	!-- Position beneath wall
849	!-- (2) - Will allways be executed.
850	!-- 'bottom and wall: use u_0,v_0 and wall functions'
851	DO i = i_left, i_right
852	DO j = j_south, j_north
853	DO k = 1, nzt
854
855	IF ( ( wall_e_x(j,i) /= 0.0_wp ).OR.( wall_e_y(j,i) /= 0.0_wp ) ) &
856	THEN
857
858	IF ( k == nzb_diff_s_inner(j,i) - 1 ) THEN
859	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
860	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
861	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
862	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
863	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
864	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
865	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
866
867	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
868	!
869	!-- Inconsistency removed: as the thermal stratification is
870	!-- not taken into account for the evaluation of the wall
871	!-- fluxes at vertical walls, the eddy viscosity km must not
872	!-- be used for the evaluation of the velocity gradients dudy
873	!-- and dwdy
874	!-- Note: The validity of the new method has not yet been
875	!-- shown, as so far no suitable data for a validation
876	!-- has been available
877	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
878	! usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
879	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
880	! wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
881	km_neutral = kappa * &
882	( usvs(k,j,i)2 + wsvs(k,j,i)2 )*0.25_wp &
883	0.5_wp * dy
884	IF ( km_neutral > 0.0_wp ) THEN
885	dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
886	dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
887	ELSE
888	dudy = 0.0_wp
889	dwdy = 0.0_wp
890	ENDIF
891	ELSE
892	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
893	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
894	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
895	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
896	ENDIF
897
898	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
899	!
900	!-- Inconsistency removed: as the thermal stratification is
901	!-- not taken into account for the evaluation of the wall
902	!-- fluxes at vertical walls, the eddy viscosity km must not
903	!-- be used for the evaluation of the velocity gradients dvdx
904	!-- and dwdx
905	!-- Note: The validity of the new method has not yet been
906	!-- shown, as so far no suitable data for a validation
907	!-- has been available
908	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
909	! vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
910	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
911	! wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
912	km_neutral = kappa * &
913	( vsus(k,j,i)2 + wsus(k,j,i)2 )*0.25_wp &
914	0.5_wp * dx
915	IF ( km_neutral > 0.0_wp ) THEN
916	dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
917	dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
918	ELSE
919	dvdx = 0.0_wp
920	dwdx = 0.0_wp
921	ENDIF
922	ELSE
923	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
924	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
925	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
926	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
927	ENDIF
928
929	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
930	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
931	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
932
933	IF ( def < 0.0_wp ) def = 0.0_wp
934
935	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
936
937	ENDIF
938	!
939	!-- (3) - will be executed only, if there is at least one level
940	!-- between (2) and (4), i.e. the topography must have a
941	!-- minimum height of 2 dz. Wall fluxes for this case have
942	!-- already been calculated for (2).
943	!-- 'wall only: use wall functions'
944
945	IF ( k >= nzb_diff_s_inner(j,i) .AND. &
946	k <= nzb_diff_s_outer(j,i)-2 ) THEN
947
948	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
949	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
950	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
951	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
952	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
953	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
954	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
955
956	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
957	!
958	!-- Inconsistency removed: as the thermal stratification
959	!-- is not taken into account for the evaluation of the
960	!-- wall fluxes at vertical walls, the eddy viscosity km
961	!-- must not be used for the evaluation of the velocity
962	!-- gradients dudy and dwdy
963	!-- Note: The validity of the new method has not yet
964	!-- been shown, as so far no suitable data for a
965	!-- validation has been available
966	km_neutral = kappa * ( usvs(k,j,i)**2 + &
967	wsvs(k,j,i)2 )0.25_wp * 0.5_wp * dy
968	IF ( km_neutral > 0.0_wp ) THEN
969	dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
970	dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
971	ELSE
972	dudy = 0.0_wp
973	dwdy = 0.0_wp
974	ENDIF
975	ELSE
976	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
977	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
978	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
979	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
980	ENDIF
981
982	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
983	!
984	!-- Inconsistency removed: as the thermal stratification
985	!-- is not taken into account for the evaluation of the
986	!-- wall fluxes at vertical walls, the eddy viscosity km
987	!-- must not be used for the evaluation of the velocity
988	!-- gradients dvdx and dwdx
989	!-- Note: The validity of the new method has not yet
990	!-- been shown, as so far no suitable data for a
991	!-- validation has been available
992	km_neutral = kappa * ( vsus(k,j,i)**2 + &
993	wsus(k,j,i)2 )0.25_wp * 0.5_wp * dx
994	IF ( km_neutral > 0.0_wp ) THEN
995	dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
996	dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
997	ELSE
998	dvdx = 0.0_wp
999	dwdx = 0.0_wp
1000	ENDIF
1001	ELSE
1002	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1003	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1004	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1005	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1006	ENDIF
1007
1008	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1009	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1010	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1011
1012	IF ( def < 0.0_wp ) def = 0.0_wp
1013
1014	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1015
1016	ENDIF
1017
1018	!
1019	!-- (4) - will allways be executed.
1020	!-- 'special case: free atmosphere' (as for case (0))
1021	IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN
1022
1023	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1024	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1025	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1026	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1027	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1028
1029	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1030	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1031	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1032	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1033	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1034
1035	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1036	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1037	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1038	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1039	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1040
1041	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1042	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1043	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1044
1045	IF ( def < 0.0_wp ) def = 0.0_wp
1046
1047	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1048
1049	ENDIF
1050
1051	ENDIF
1052
1053	ENDDO
1054	ENDDO
1055	ENDDO
1056
1057	!
1058	!-- Position without adjacent wall
1059	!-- (1) - will allways be executed.
1060	!-- 'bottom only: use u_0,v_0'
1061	DO i = i_left, i_right
1062	DO j = j_south, j_north
1063	DO k = 1, nzt
1064
1065	IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
1066	THEN
1067
1068	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
1069
1070	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1071	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1072	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1073	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1074	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
1075
1076	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1077	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1078	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1079	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1080	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
1081
1082	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1083	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1084	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1085	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1086	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1087
1088	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1089	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1090	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1091
1092	IF ( def < 0.0_wp ) def = 0.0_wp
1093
1094	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1095
1096	ENDIF
1097
1098	ENDIF
1099
1100	ENDDO
1101	ENDDO
1102	ENDDO
1103
1104	ELSEIF ( use_surface_fluxes ) THEN
1105
1106	DO i = i_left, i_right
1107	DO j = j_south, j_north
1108	DO k = 1, nzt
1109
1110	IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN
1111
1112	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1113	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1114	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1115	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1116	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1117
1118	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1119	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1120	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1121	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1122	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1123
1124	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1125	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1126	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1127	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1128	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1129
1130	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1131	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1132	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1133
1134	IF ( def < 0.0_wp ) def = 0.0_wp
1135
1136	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1137
1138	ENDIF
1139
1140	ENDDO
1141	ENDDO
1142	ENDDO
1143
1144	ENDIF
1145
1146	!
1147	!-- If required, calculate TKE production by buoyancy
1148	IF ( .NOT. neutral ) THEN
1149
1150	IF ( .NOT. humidity ) THEN
1151
1152	IF ( use_single_reference_value ) THEN
1153
1154	IF ( ocean ) THEN
1155	!
1156	!-- So far in the ocean no special treatment of density flux
1157	!-- in the bottom and top surface layer
1158	DO i = i_left, i_right
1159	DO j = j_south, j_north
1160	DO k = 1, nzt
1161	IF ( k > nzb_s_inner(j,i) ) THEN
1162	tend(k,j,i) = tend(k,j,i) + &
1163	kh(k,j,i) * g / rho_reference * &
1164	( rho(k+1,j,i) - rho(k-1,j,i) ) * &
1165	dd2zu(k)
1166	ENDIF
1167	ENDDO
1168	ENDDO
1169	ENDDO
1170
1171	ELSE
1172
1173	DO i = i_left, i_right
1174	DO j = j_south, j_north
1175	DO k = 1, nzt_diff
1176	IF ( k >= nzb_diff_s_inner(j,i) ) THEN
1177	tend(k,j,i) = tend(k,j,i) - &
1178	kh(k,j,i) * g / pt_reference * &
1179	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
1180	dd2zu(k)
1181	ENDIF
1182
1183	IF ( k == nzb_diff_s_inner(j,i)-1 .AND. &
1184	use_surface_fluxes ) THEN
1185	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
1186	shf(j,i)
1187	ENDIF
1188
1189	IF ( k == nzt .AND. use_top_fluxes ) THEN
1190	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
1191	tswst(j,i)
1192	ENDIF
1193	ENDDO
1194	ENDDO
1195	ENDDO
1196
1197	ENDIF
1198
1199	ELSE
1200
1201	IF ( ocean ) THEN
1202	!
1203	!-- So far in the ocean no special treatment of density flux
1204	!-- in the bottom and top surface layer
1205	DO i = i_left, i_right
1206	DO j = j_south, j_north
1207	DO k = 1, nzt
1208	IF ( k > nzb_s_inner(j,i) ) THEN
1209	tend(k,j,i) = tend(k,j,i) + &
1210	kh(k,j,i) * g / rho(k,j,i) * &
1211	( rho(k+1,j,i) - rho(k-1,j,i) ) * &
1212	dd2zu(k)
1213	ENDIF
1214	ENDDO
1215	ENDDO
1216	ENDDO
1217
1218	ELSE
1219
1220	DO i = i_left, i_right
1221	DO j = j_south, j_north
1222	DO k = 1, nzt_diff
1223	IF( k >= nzb_diff_s_inner(j,i) ) THEN
1224	tend(k,j,i) = tend(k,j,i) - &
1225	kh(k,j,i) * g / pt(k,j,i) * &
1226	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
1227	dd2zu(k)
1228	ENDIF
1229
1230	IF ( k == nzb_diff_s_inner(j,i)-1 .AND. &
1231	use_surface_fluxes ) THEN
1232	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
1233	shf(j,i)
1234	ENDIF
1235
1236	IF ( k == nzt .AND. use_top_fluxes ) THEN
1237	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
1238	tswst(j,i)
1239	ENDIF
1240	ENDDO
1241	ENDDO
1242	ENDDO
1243
1244	ENDIF
1245
1246	ENDIF
1247
1248	ELSE
1249	!
1250	!++ This part gives the PGI compiler problems in the previous loop
1251	!++ even without any acc statements????
1252	! STOP '+++ production_e problems with acc-directives'
1253	! !acc loop
1254	! DO i = nxl, nxr
1255	! DO j = nys, nyn
1256	! !acc loop vector
1257	! DO k = 1, nzt_diff
1258	!
1259	! IF ( k >= nzb_diff_s_inner(j,i) ) THEN
1260	!
1261	! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1262	! k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1263	! k2 = 0.61_wp * pt(k,j,i)
1264	! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
1265	! g / vpt(k,j,i) * &
1266	! ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1267	! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
1268	! ) * dd2zu(k)
1269	! ELSE IF ( cloud_physics ) THEN
1270	! IF ( ql(k,j,i) == 0.0_wp ) THEN
1271	! k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1272	! k2 = 0.61_wp * pt(k,j,i)
1273	! ELSE
1274	! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1275	! temp = theta * t_d_pt(k)
1276	! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1277	! ( q(k,j,i) - ql(k,j,i) ) * &
1278	! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1279	! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1280	! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1281	! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1282	! ENDIF
1283	! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
1284	! g / vpt(k,j,i) * &
1285	! ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1286	! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
1287	! ) * dd2zu(k)
1288	! ELSE IF ( cloud_droplets ) THEN
1289	! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1290	! k2 = 0.61_wp * pt(k,j,i)
1291	! tend(k,j,i) = tend(k,j,i) - &
1292	! kh(k,j,i) * g / vpt(k,j,i) * &
1293	! ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &
1294	! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
1295	! pt(k,j,i) * ( ql(k+1,j,i) - &
1296	! ql(k-1,j,i) ) ) * dd2zu(k)
1297	! ENDIF
1298	!
1299	! ENDIF
1300	!
1301	! ENDDO
1302	! ENDDO
1303	! ENDDO
1304	!
1305
1306	!!++ Next two loops are probably very inefficiently parallellized
1307	!!++ and will require better optimization
1308	! IF ( use_surface_fluxes ) THEN
1309	!
1310	! !acc loop
1311	! DO i = nxl, nxr
1312	! DO j = nys, nyn
1313	! !acc loop vector
1314	! DO k = 1, nzt_diff
1315	!
1316	! IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
1317	!
1318	! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1319	! k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1320	! k2 = 0.61_wp * pt(k,j,i)
1321	! ELSE IF ( cloud_physics ) THEN
1322	! IF ( ql(k,j,i) == 0.0_wp ) THEN
1323	! k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1324	! k2 = 0.61_wp * pt(k,j,i)
1325	! ELSE
1326	! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1327	! temp = theta * t_d_pt(k)
1328	! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1329	! ( q(k,j,i) - ql(k,j,i) ) * &
1330	! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&
1331	! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1332	! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1333	! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1334	! ENDIF
1335	! ELSE IF ( cloud_droplets ) THEN
1336	! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1337	! k2 = 0.61_wp * pt(k,j,i)
1338	! ENDIF
1339	!
1340	! tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1341	! ( k1* shf(j,i) + k2 * qsws(j,i) )
1342	! ENDIF
1343	!
1344	! ENDDO
1345	! ENDDO
1346	! ENDDO
1347	!
1348	! ENDIF
1349	!
1350	! IF ( use_top_fluxes ) THEN
1351	!
1352	! !acc loop
1353	! DO i = nxl, nxr
1354	! DO j = nys, nyn
1355	! !acc loop vector
1356	! DO k = 1, nzt
1357	! IF ( k == nzt ) THEN
1358	!
1359	! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1360	! k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1361	! k2 = 0.61_wp * pt(k,j,i)
1362	! ELSE IF ( cloud_physics ) THEN
1363	! IF ( ql(k,j,i) == 0.0_wp ) THEN
1364	! k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1365	! k2 = 0.61_wp * pt(k,j,i)
1366	! ELSE
1367	! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1368	! temp = theta * t_d_pt(k)
1369	! k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1370	! ( q(k,j,i) - ql(k,j,i) ) * &
1371	! ( 1.0_wp + 0.622_wp * l_d_r / temp ) ) /&
1372	! ( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1373	! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1374	! k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1375	! ENDIF
1376	! ELSE IF ( cloud_droplets ) THEN
1377	! k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1378	! k2 = 0.61_wp * pt(k,j,i)
1379	! ENDIF
1380	!
1381	! tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1382	! ( k1* tswst(j,i) + k2 * qswst(j,i) )
1383	!
1384	! ENDIF
1385	!
1386	! ENDDO
1387	! ENDDO
1388	! ENDDO
1389	!
1390	! ENDIF
1391
1392	ENDIF
1393
1394	ENDIF
1395	!$acc end kernels
1396
1397	END SUBROUTINE production_e_acc
1398
1399
1400	!------------------------------------------------------------------------------!
1401	! Description:
1402	! ------------
1403	!> Call for grid point i,j
1404	!------------------------------------------------------------------------------!
1405	SUBROUTINE production_e_ij( i, j )
1406
1407	USE arrays_3d, &
1408	ONLY: ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho, shf, &
1409	tend, tswst, u, v, vpt, w
1410
1411	USE cloud_parameters, &
1412	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
1413
1414	USE control_parameters, &
1415	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
1416	humidity, kappa, neutral, ocean, pt_reference, &
1417	rho_reference, use_single_reference_value, &
1418	use_surface_fluxes, use_top_fluxes
1419
1420	USE grid_variables, &
1421	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
1422
1423	USE indices, &
1424	ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, &
1425	nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
1426
1427	IMPLICIT NONE
1428
1429	INTEGER(iwp) :: i !<
1430	INTEGER(iwp) :: j !<
1431	INTEGER(iwp) :: k !<
1432
1433	REAL(wp) :: def !<
1434	REAL(wp) :: dudx !<
1435	REAL(wp) :: dudy !<
1436	REAL(wp) :: dudz !<
1437	REAL(wp) :: dvdx !<
1438	REAL(wp) :: dvdy !<
1439	REAL(wp) :: dvdz !<
1440	REAL(wp) :: dwdx !<
1441	REAL(wp) :: dwdy !<
1442	REAL(wp) :: dwdz !<
1443	REAL(wp) :: k1 !<
1444	REAL(wp) :: k2 !<
1445	REAL(wp) :: km_neutral !<
1446	REAL(wp) :: theta !<
1447	REAL(wp) :: temp !<
1448
1449	REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !<
1450	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
1451	REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !<
1452	REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !<
1453
1454	!
1455	!-- Calculate TKE production by shear
1456	DO k = nzb_diff_s_outer(j,i), nzt
1457
1458	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1459	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1460	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1461	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1462	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1463
1464	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1465	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1466	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1467	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1468	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1469
1470	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1471	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1472	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1473	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1474	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1475
1476	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) &
1477	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
1478	+ 2.0_wp * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1479
1480	IF ( def < 0.0_wp ) def = 0.0_wp
1481
1482	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1483
1484	ENDDO
1485
1486	IF ( constant_flux_layer ) THEN
1487
1488	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) THEN
1489
1490	!
1491	!-- Position beneath wall
1492	!-- (2) - Will allways be executed.
1493	!-- 'bottom and wall: use u_0,v_0 and wall functions'
1494	k = nzb_diff_s_inner(j,i)-1
1495
1496	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1497	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1498	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
1499	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1500	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1501	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
1502	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1503
1504	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
1505	!
1506	!-- Inconsistency removed: as the thermal stratification
1507	!-- is not taken into account for the evaluation of the
1508	!-- wall fluxes at vertical walls, the eddy viscosity km
1509	!-- must not be used for the evaluation of the velocity
1510	!-- gradients dudy and dwdy
1511	!-- Note: The validity of the new method has not yet
1512	!-- been shown, as so far no suitable data for a
1513	!-- validation has been available
1514	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
1515	usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
1516	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
1517	wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
1518	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25_wp &
1519	0.5_wp * dy
1520	IF ( km_neutral > 0.0_wp ) THEN
1521	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
1522	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
1523	ELSE
1524	dudy = 0.0_wp
1525	dwdy = 0.0_wp
1526	ENDIF
1527	ELSE
1528	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1529	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1530	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1531	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1532	ENDIF
1533
1534	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
1535	!
1536	!-- Inconsistency removed: as the thermal stratification
1537	!-- is not taken into account for the evaluation of the
1538	!-- wall fluxes at vertical walls, the eddy viscosity km
1539	!-- must not be used for the evaluation of the velocity
1540	!-- gradients dvdx and dwdx
1541	!-- Note: The validity of the new method has not yet
1542	!-- been shown, as so far no suitable data for a
1543	!-- validation has been available
1544	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
1545	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
1546	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
1547	wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
1548	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25_wp &
1549	0.5_wp * dx
1550	IF ( km_neutral > 0.0_wp ) THEN
1551	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
1552	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
1553	ELSE
1554	dvdx = 0.0_wp
1555	dwdx = 0.0_wp
1556	ENDIF
1557	ELSE
1558	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1559	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1560	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1561	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1562	ENDIF
1563
1564	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1565	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1566	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1567
1568	IF ( def < 0.0_wp ) def = 0.0_wp
1569
1570	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1571
1572	!
1573	!-- (3) - will be executed only, if there is at least one level
1574	!-- between (2) and (4), i.e. the topography must have a
1575	!-- minimum height of 2 dz. Wall fluxes for this case have
1576	!-- already been calculated for (2).
1577	!-- 'wall only: use wall functions'
1578	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
1579
1580	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1581	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1582	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1583	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1584	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1585	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1586	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1587
1588	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
1589	!
1590	!-- Inconsistency removed: as the thermal stratification
1591	!-- is not taken into account for the evaluation of the
1592	!-- wall fluxes at vertical walls, the eddy viscosity km
1593	!-- must not be used for the evaluation of the velocity
1594	!-- gradients dudy and dwdy
1595	!-- Note: The validity of the new method has not yet
1596	!-- been shown, as so far no suitable data for a
1597	!-- validation has been available
1598	km_neutral = kappa * ( usvs(k)**2 + &
1599	wsvs(k)2 )0.25_wp * 0.5_wp * dy
1600	IF ( km_neutral > 0.0_wp ) THEN
1601	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
1602	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
1603	ELSE
1604	dudy = 0.0_wp
1605	dwdy = 0.0_wp
1606	ENDIF
1607	ELSE
1608	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1609	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1610	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1611	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1612	ENDIF
1613
1614	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
1615	!
1616	!-- Inconsistency removed: as the thermal stratification
1617	!-- is not taken into account for the evaluation of the
1618	!-- wall fluxes at vertical walls, the eddy viscosity km
1619	!-- must not be used for the evaluation of the velocity
1620	!-- gradients dvdx and dwdx
1621	!-- Note: The validity of the new method has not yet
1622	!-- been shown, as so far no suitable data for a
1623	!-- validation has been available
1624	km_neutral = kappa * ( vsus(k)**2 + &
1625	wsus(k)2 )0.25_wp * 0.5_wp * dx
1626	IF ( km_neutral > 0.0_wp ) THEN
1627	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
1628	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
1629	ELSE
1630	dvdx = 0.0_wp
1631	dwdx = 0.0_wp
1632	ENDIF
1633	ELSE
1634	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1635	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1636	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1637	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1638	ENDIF
1639
1640	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1641	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1642	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1643
1644	IF ( def < 0.0_wp ) def = 0.0_wp
1645
1646	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1647
1648	ENDDO
1649
1650	!
1651	!-- (4) - will allways be executed.
1652	!-- 'special case: free atmosphere' (as for case (0))
1653	k = nzb_diff_s_outer(j,i)-1
1654
1655	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1656	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1657	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1658	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1659	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1660
1661	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1662	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1663	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1664	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1665	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1666
1667	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1668	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1669	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1670	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1671	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1672
1673	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1674	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1675	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1676
1677	IF ( def < 0.0_wp ) def = 0.0_wp
1678
1679	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1680
1681	ELSE
1682
1683	!
1684	!-- Position without adjacent wall
1685	!-- (1) - will allways be executed.
1686	!-- 'bottom only: use u_0,v_0'
1687	k = nzb_diff_s_inner(j,i)-1
1688
1689	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1690	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1691	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1692	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1693	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
1694
1695	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1696	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1697	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1698	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1699	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
1700
1701	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1702	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1703	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1704	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1705	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1706
1707	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) &
1708	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
1709	+ 2.0_wp * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1710
1711	IF ( def < 0.0_wp ) def = 0.0_wp
1712
1713	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1714
1715	ENDIF
1716
1717	ELSEIF ( use_surface_fluxes ) THEN
1718
1719	k = nzb_diff_s_outer(j,i)-1
1720
1721	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1722	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1723	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1724	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1725	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1726
1727	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1728	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1729	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1730	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1731	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1732
1733	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1734	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1735	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1736	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1737	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1738
1739	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1740	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1741	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1742
1743	IF ( def < 0.0_wp ) def = 0.0_wp
1744
1745	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1746
1747	ENDIF
1748
1749	!
1750	!-- If required, calculate TKE production by buoyancy
1751	IF ( .NOT. neutral ) THEN
1752
1753	IF ( .NOT. humidity ) THEN
1754
1755	IF ( use_single_reference_value ) THEN
1756
1757	IF ( ocean ) THEN
1758	!
1759	!-- So far in the ocean no special treatment of density flux in
1760	!-- the bottom and top surface layer
1761	DO k = nzb_s_inner(j,i)+1, nzt
1762	tend(k,j,i) = tend(k,j,i) + &
1763	kh(k,j,i) * g / rho_reference * &
1764	( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
1765	ENDDO
1766
1767	ELSE
1768
1769	DO k = nzb_diff_s_inner(j,i), nzt_diff
1770	tend(k,j,i) = tend(k,j,i) - &
1771	kh(k,j,i) * g / pt_reference * &
1772	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
1773	ENDDO
1774
1775	IF ( use_surface_fluxes ) THEN
1776	k = nzb_diff_s_inner(j,i)-1
1777	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
1778	ENDIF
1779
1780	IF ( use_top_fluxes ) THEN
1781	k = nzt
1782	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
1783	ENDIF
1784
1785	ENDIF
1786
1787	ELSE
1788
1789	IF ( ocean ) THEN
1790	!
1791	!-- So far in the ocean no special treatment of density flux in
1792	!-- the bottom and top surface layer
1793	DO k = nzb_s_inner(j,i)+1, nzt
1794	tend(k,j,i) = tend(k,j,i) + &
1795	kh(k,j,i) * g / rho(k,j,i) * &
1796	( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
1797	ENDDO
1798
1799	ELSE
1800
1801	DO k = nzb_diff_s_inner(j,i), nzt_diff
1802	tend(k,j,i) = tend(k,j,i) - &
1803	kh(k,j,i) * g / pt(k,j,i) * &
1804	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
1805	ENDDO
1806
1807	IF ( use_surface_fluxes ) THEN
1808	k = nzb_diff_s_inner(j,i)-1
1809	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
1810	ENDIF
1811
1812	IF ( use_top_fluxes ) THEN
1813	k = nzt
1814	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
1815	ENDIF
1816
1817	ENDIF
1818
1819	ENDIF
1820
1821	ELSE
1822
1823	DO k = nzb_diff_s_inner(j,i), nzt_diff
1824
1825	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1826	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1827	k2 = 0.61_wp * pt(k,j,i)
1828	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
1829	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1830	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
1831	) * dd2zu(k)
1832	ELSE IF ( cloud_physics ) THEN
1833	IF ( ql(k,j,i) == 0.0_wp ) THEN
1834	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1835	k2 = 0.61_wp * pt(k,j,i)
1836	ELSE
1837	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1838	temp = theta * t_d_pt(k)
1839	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1840	( q(k,j,i) - ql(k,j,i) ) * &
1841	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1842	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1843	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1844	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1845	ENDIF
1846	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
1847	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1848	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
1849	) * dd2zu(k)
1850	ELSE IF ( cloud_droplets ) THEN
1851	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1852	k2 = 0.61_wp * pt(k,j,i)
1853	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
1854	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1855	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
1856	pt(k,j,i) * ( ql(k+1,j,i) - &
1857	ql(k-1,j,i) ) ) * dd2zu(k)
1858	ENDIF
1859	ENDDO
1860
1861	IF ( use_surface_fluxes ) THEN
1862	k = nzb_diff_s_inner(j,i)-1
1863
1864	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1865	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1866	k2 = 0.61_wp * pt(k,j,i)
1867	ELSE IF ( cloud_physics ) THEN
1868	IF ( ql(k,j,i) == 0.0_wp ) THEN
1869	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1870	k2 = 0.61_wp * pt(k,j,i)
1871	ELSE
1872	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1873	temp = theta * t_d_pt(k)
1874	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1875	( q(k,j,i) - ql(k,j,i) ) * &
1876	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1877	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1878	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1879	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1880	ENDIF
1881	ELSE IF ( cloud_droplets ) THEN
1882	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1883	k2 = 0.61_wp * pt(k,j,i)
1884	ENDIF
1885
1886	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1887	( k1* shf(j,i) + k2 * qsws(j,i) )
1888	ENDIF
1889
1890	IF ( use_top_fluxes ) THEN
1891	k = nzt
1892
1893	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1894	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1895	k2 = 0.61_wp * pt(k,j,i)
1896	ELSE IF ( cloud_physics ) THEN
1897	IF ( ql(k,j,i) == 0.0_wp ) THEN
1898	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1899	k2 = 0.61_wp * pt(k,j,i)
1900	ELSE
1901	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1902	temp = theta * t_d_pt(k)
1903	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1904	( q(k,j,i) - ql(k,j,i) ) * &
1905	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1906	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1907	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1908	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1909	ENDIF
1910	ELSE IF ( cloud_droplets ) THEN
1911	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1912	k2 = 0.61_wp * pt(k,j,i)
1913	ENDIF
1914
1915	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1916	( k1* tswst(j,i) + k2 * qswst(j,i) )
1917	ENDIF
1918
1919	ENDIF
1920
1921	ENDIF
1922
1923	END SUBROUTINE production_e_ij
1924
1925
1926	!------------------------------------------------------------------------------!
1927	! Description:
1928	! ------------
1929	!> @todo Missing subroutine description.
1930	!------------------------------------------------------------------------------!
1931	SUBROUTINE production_e_init
1932
1933	USE arrays_3d, &
1934	ONLY: kh, km, u, us, usws, v, vsws, zu
1935
1936	USE control_parameters, &
1937	ONLY: constant_flux_layer, kappa
1938
1939	USE indices, &
1940	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner, &
1941	nzb_v_inner
1942
1943	IMPLICIT NONE
1944
1945	INTEGER(iwp) :: i !<
1946	INTEGER(iwp) :: j !<
1947	INTEGER(iwp) :: ku !<
1948	INTEGER(iwp) :: kv !<
1949
1950	IF ( constant_flux_layer ) THEN
1951
1952	IF ( first_call ) THEN
1953	ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
1954	u_0 = 0.0_wp ! just to avoid access of uninitialized memory
1955	v_0 = 0.0_wp ! within exchange_horiz_2d
1956	first_call = .FALSE.
1957	ENDIF
1958
1959	!
1960	!-- Calculate a virtual velocity at the surface in a way that the
1961	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1962	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1963	!-- production term at k=1 (see production_e_ij).
1964	!-- The velocity gradient has to be limited in case of too small km
1965	!-- (otherwise the timestep may be significantly reduced by large
1966	!-- surface winds).
1967	!-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are
1968	!-- not available in case of non-cyclic boundary conditions.
1969	!-- WARNING: the exact analytical solution would require the determination
1970	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
1971	!$OMP PARALLEL DO PRIVATE( ku, kv )
1972	DO i = nxl, nxr+1
1973	DO j = nys, nyn+1
1974
1975	ku = nzb_u_inner(j,i)+1
1976	kv = nzb_v_inner(j,i)+1
1977
1978	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
1979	( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) + &
1980	1.0E-20_wp )
1981	! ( us(j,i) * kappa * zu(1) )
1982	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
1983	( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) + &
1984	1.0E-20_wp )
1985	! ( us(j,i) * kappa * zu(1) )
1986
1987	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
1988	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
1989	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
1990	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
1991
1992	ENDDO
1993	ENDDO
1994
1995	CALL exchange_horiz_2d( u_0 )
1996	CALL exchange_horiz_2d( v_0 )
1997
1998	ENDIF
1999
2000	END SUBROUTINE production_e_init
2001
2002	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |