Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |