Home

Context Navigation

production_e.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 10 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:

openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:

bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

Property svn:keywords set to Id

File size: 82.7 KB

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

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

Context Navigation

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

Download in other formats: