Home

Context Navigation

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

Last change on this file since 2193 was 2127, checked in by raasch, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 57.9 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-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: production_e.f90 2127 2017-01-20 16:06:40Z raasch $
27	!
28	! 2126 2017-01-20 15:54:21Z raasch
29	! density in ocean case replaced by potential density
30	!
31	! 2118 2017-01-17 16:38:49Z raasch
32	! OpenACC version of subroutine removed
33	!
34	! 2031 2016-10-21 15:11:58Z knoop
35	! renamed variable rho to rho_ocean
36	!
37	! 2000 2016-08-20 18:09:15Z knoop
38	! Forced header and separation lines into 80 columns
39	!
40	! 1873 2016-04-18 14:50:06Z maronga
41	! Module renamed (removed _mod)
42	!
43	!
44	! 1850 2016-04-08 13:29:27Z maronga
45	! Module renamed
46	!
47	!
48	! 1691 2015-10-26 16:17:44Z maronga
49	! Renamed prandtl_layer to constant_flux_layer.
50	!
51	! 1682 2015-10-07 23:56:08Z knoop
52	! Code annotations made doxygen readable
53	!
54	! 1374 2014-04-25 12:55:07Z raasch
55	! nzb_s_outer removed from acc-present-list
56	!
57	! 1353 2014-04-08 15:21:23Z heinze
58	! REAL constants provided with KIND-attribute
59	!
60	! 1342 2014-03-26 17:04:47Z kanani
61	! REAL constants defined as wp-kind
62	!
63	! 1320 2014-03-20 08:40:49Z raasch
64	! ONLY-attribute added to USE-statements,
65	! kind-parameters added to all INTEGER and REAL declaration statements,
66	! kinds are defined in new module kinds,
67	! old module precision_kind is removed,
68	! revision history before 2012 removed,
69	! comment fields (!:) to be used for variable explanations added to
70	! all variable declaration statements
71	!
72	! 1257 2013-11-08 15:18:40Z raasch
73	! openacc loop and loop vector clauses removed, declare create moved after
74	! the FORTRAN declaration statement
75	!
76	! 1179 2013-06-14 05:57:58Z raasch
77	! use_reference renamed use_single_reference_value
78	!
79	! 1128 2013-04-12 06:19:32Z raasch
80	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
81	! j_north
82	!
83	! 1036 2012-10-22 13:43:42Z raasch
84	! code put under GPL (PALM 3.9)
85	!
86	! 1015 2012-09-27 09:23:24Z raasch
87	! accelerator version (*_acc) added
88	!
89	! 1007 2012-09-19 14:30:36Z franke
90	! Bugfix: calculation of buoyancy production has to consider the liquid water
91	! mixing ratio in case of cloud droplets
92	!
93	! 940 2012-07-09 14:31:00Z raasch
94	! TKE production by buoyancy can be switched off in case of runs with pure
95	! neutral stratification
96	!
97	! Revision 1.1 1997/09/19 07:45:35 raasch
98	! Initial revision
99	!
100	!
101	! Description:
102	! ------------
103	!> Production terms (shear + buoyancy) of the TKE.
104	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
105	!> not considered well!
106	!------------------------------------------------------------------------------!
107	MODULE production_e_mod
108
109
110	USE wall_fluxes_mod, &
111	ONLY: wall_fluxes_e
112
113	USE kinds
114
115	PRIVATE
116	PUBLIC production_e, production_e_init
117
118	LOGICAL, SAVE :: first_call = .TRUE. !<
119
120	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0 !<
121	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: v_0 !<
122
123	INTERFACE production_e
124	MODULE PROCEDURE production_e
125	MODULE PROCEDURE production_e_ij
126	END INTERFACE production_e
127
128	INTERFACE production_e_init
129	MODULE PROCEDURE production_e_init
130	END INTERFACE production_e_init
131
132	CONTAINS
133
134
135	!------------------------------------------------------------------------------!
136	! Description:
137	! ------------
138	!> Call for all grid points
139	!------------------------------------------------------------------------------!
140	SUBROUTINE production_e
141
142	USE arrays_3d, &
143	ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf, &
144	tend, tswst, u, v, vpt, w
145
146	USE cloud_parameters, &
147	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
148
149	USE control_parameters, &
150	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
151	humidity, kappa, neutral, ocean, pt_reference, &
152	rho_reference, use_single_reference_value, &
153	use_surface_fluxes, use_top_fluxes
154
155	USE grid_variables, &
156	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
157
158	USE indices, &
159	ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, &
160	nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
161
162	IMPLICIT NONE
163
164	INTEGER(iwp) :: i !<
165	INTEGER(iwp) :: j !<
166	INTEGER(iwp) :: k !<
167
168	REAL(wp) :: def !<
169	REAL(wp) :: dudx !<
170	REAL(wp) :: dudy !<
171	REAL(wp) :: dudz !<
172	REAL(wp) :: dvdx !<
173	REAL(wp) :: dvdy !<
174	REAL(wp) :: dvdz !<
175	REAL(wp) :: dwdx !<
176	REAL(wp) :: dwdy !<
177	REAL(wp) :: dwdz !<
178	REAL(wp) :: k1 !<
179	REAL(wp) :: k2 !<
180	REAL(wp) :: km_neutral !<
181	REAL(wp) :: theta !<
182	REAL(wp) :: temp !<
183
184	! REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs
185	REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !<
186	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
187	REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !<
188	REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !<
189
190	!
191	!-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
192	!-- vertical walls, if neccessary
193	!-- So far, results are slightly different from the ij-Version.
194	!-- Therefore, ij-Version is called further below within the ij-loops.
195	! IF ( topography /= 'flat' ) THEN
196	! CALL wall_fluxes_e( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )
197	! CALL wall_fluxes_e( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )
198	! CALL wall_fluxes_e( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )
199	! CALL wall_fluxes_e( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )
200	! ENDIF
201
202
203	DO i = nxl, nxr
204
205	!
206	!-- Calculate TKE production by shear
207	DO j = nys, nyn
208	DO k = nzb_diff_s_outer(j,i), nzt
209
210	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
211	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
212	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
213	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
214	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
215
216	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
217	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
218	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
219	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
220	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
221
222	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
223	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
224	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
225	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
226	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
227
228	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
229	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
230	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
231
232	IF ( def < 0.0_wp ) def = 0.0_wp
233
234	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
235
236	ENDDO
237	ENDDO
238
239	IF ( constant_flux_layer ) THEN
240
241	!
242	!-- Position beneath wall
243	!-- (2) - Will allways be executed.
244	!-- 'bottom and wall: use u_0,v_0 and wall functions'
245	DO j = nys, nyn
246
247	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
248	THEN
249
250	k = nzb_diff_s_inner(j,i) - 1
251	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
252	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
253	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
254	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
255	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
256	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
257	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
258
259	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
260	!
261	!-- Inconsistency removed: as the thermal stratification is
262	!-- not taken into account for the evaluation of the wall
263	!-- fluxes at vertical walls, the eddy viscosity km must not
264	!-- be used for the evaluation of the velocity gradients dudy
265	!-- and dwdy
266	!-- Note: The validity of the new method has not yet been
267	!-- shown, as so far no suitable data for a validation
268	!-- has been available
269	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
270	usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
271	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
272	wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
273	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25_wp &
274	0.5_wp * dy
275	IF ( km_neutral > 0.0_wp ) THEN
276	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
277	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
278	ELSE
279	dudy = 0.0_wp
280	dwdy = 0.0_wp
281	ENDIF
282	ELSE
283	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
284	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
285	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
286	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
287	ENDIF
288
289	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
290	!
291	!-- Inconsistency removed: as the thermal stratification is
292	!-- not taken into account for the evaluation of the wall
293	!-- fluxes at vertical walls, the eddy viscosity km must not
294	!-- be used for the evaluation of the velocity gradients dvdx
295	!-- and dwdx
296	!-- Note: The validity of the new method has not yet been
297	!-- shown, as so far no suitable data for a validation
298	!-- has been available
299	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
300	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
301	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
302	wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
303	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25_wp &
304	0.5_wp * dx
305	IF ( km_neutral > 0.0_wp ) THEN
306	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
307	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
308	ELSE
309	dvdx = 0.0_wp
310	dwdx = 0.0_wp
311	ENDIF
312	ELSE
313	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
314	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
315	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
316	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
317	ENDIF
318
319	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
320	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
321	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
322
323	IF ( def < 0.0_wp ) def = 0.0_wp
324
325	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
326
327
328	!
329	!-- (3) - will be executed only, if there is at least one level
330	!-- between (2) and (4), i.e. the topography must have a
331	!-- minimum height of 2 dz. Wall fluxes for this case have
332	!-- already been calculated for (2).
333	!-- 'wall only: use wall functions'
334
335	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
336
337	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
338	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
339	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
340	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
341	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
342	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
343	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
344
345	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
346	!
347	!-- Inconsistency removed: as the thermal stratification
348	!-- is not taken into account for the evaluation of the
349	!-- wall fluxes at vertical walls, the eddy viscosity km
350	!-- must not be used for the evaluation of the velocity
351	!-- gradients dudy and dwdy
352	!-- Note: The validity of the new method has not yet
353	!-- been shown, as so far no suitable data for a
354	!-- validation has been available
355	km_neutral = kappa * ( usvs(k)**2 + &
356	wsvs(k)2 )0.25_wp * 0.5_wp * dy
357	IF ( km_neutral > 0.0_wp ) THEN
358	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
359	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
360	ELSE
361	dudy = 0.0_wp
362	dwdy = 0.0_wp
363	ENDIF
364	ELSE
365	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
366	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
367	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
368	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
369	ENDIF
370
371	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
372	!
373	!-- Inconsistency removed: as the thermal stratification
374	!-- is not taken into account for the evaluation of the
375	!-- wall fluxes at vertical walls, the eddy viscosity km
376	!-- must not be used for the evaluation of the velocity
377	!-- gradients dvdx and dwdx
378	!-- Note: The validity of the new method has not yet
379	!-- been shown, as so far no suitable data for a
380	!-- validation has been available
381	km_neutral = kappa * ( vsus(k)**2 + &
382	wsus(k)2 )0.25_wp * 0.5_wp * dx
383	IF ( km_neutral > 0.0_wp ) THEN
384	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
385	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
386	ELSE
387	dvdx = 0.0_wp
388	dwdx = 0.0_wp
389	ENDIF
390	ELSE
391	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
392	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
393	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
394	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
395	ENDIF
396
397	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
398	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
399	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
400
401	IF ( def < 0.0_wp ) def = 0.0_wp
402
403	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
404
405	ENDDO
406
407	ENDIF
408
409	ENDDO
410
411	!
412	!-- (4) - will allways be executed.
413	!-- 'special case: free atmosphere' (as for case (0))
414	DO j = nys, nyn
415
416	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
417	THEN
418
419	k = nzb_diff_s_outer(j,i)-1
420
421	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
422	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
423	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
424	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
425	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
426
427	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
428	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
429	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
430	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
431	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
432
433	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
434	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
435	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
436	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
437	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
438
439	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
440	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
441	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
442
443	IF ( def < 0.0_wp ) def = 0.0_wp
444
445	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
446
447	ENDIF
448
449	ENDDO
450
451	!
452	!-- Position without adjacent wall
453	!-- (1) - will allways be executed.
454	!-- 'bottom only: use u_0,v_0'
455	DO j = nys, nyn
456
457	IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
458	THEN
459
460	k = nzb_diff_s_inner(j,i)-1
461
462	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
463	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
464	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
465	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
466	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
467
468	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
469	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
470	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
471	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
472	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
473
474	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
475	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
476	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
477	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
478	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
479
480	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
481	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
482	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
483
484	IF ( def < 0.0_wp ) def = 0.0_wp
485
486	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
487
488	ENDIF
489
490	ENDDO
491
492	ELSEIF ( use_surface_fluxes ) THEN
493
494	DO j = nys, nyn
495
496	k = nzb_diff_s_outer(j,i)-1
497
498	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
499	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
500	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
501	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
502	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
503
504	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
505	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
506	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
507	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
508	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
509
510	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
511	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
512	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
513	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
514	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
515
516	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
517	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
518	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
519
520	IF ( def < 0.0_wp ) def = 0.0_wp
521
522	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
523
524	ENDDO
525
526	ENDIF
527
528	!
529	!-- If required, calculate TKE production by buoyancy
530	IF ( .NOT. neutral ) THEN
531
532	IF ( .NOT. humidity ) THEN
533
534	IF ( use_single_reference_value ) THEN
535
536	IF ( ocean ) THEN
537	!
538	!-- So far in the ocean no special treatment of density flux
539	!-- in the bottom and top surface layer
540	DO j = nys, nyn
541	DO k = nzb_s_inner(j,i)+1, nzt
542	tend(k,j,i) = tend(k,j,i) + &
543	kh(k,j,i) * g / rho_reference * &
544	( prho(k+1,j,i) - prho(k-1,j,i) ) * &
545	dd2zu(k)
546	ENDDO
547	ENDDO
548
549	ELSE
550
551	DO j = nys, nyn
552	DO k = nzb_diff_s_inner(j,i), nzt_diff
553	tend(k,j,i) = tend(k,j,i) - &
554	kh(k,j,i) * g / pt_reference * &
555	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
556	dd2zu(k)
557	ENDDO
558
559	IF ( use_surface_fluxes ) THEN
560	k = nzb_diff_s_inner(j,i)-1
561	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
562	shf(j,i)
563	ENDIF
564
565	IF ( use_top_fluxes ) THEN
566	k = nzt
567	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
568	tswst(j,i)
569	ENDIF
570	ENDDO
571
572	ENDIF
573
574	ELSE
575
576	IF ( ocean ) THEN
577	!
578	!-- So far in the ocean no special treatment of density flux
579	!-- in the bottom and top surface layer
580	DO j = nys, nyn
581	DO k = nzb_s_inner(j,i)+1, nzt
582	tend(k,j,i) = tend(k,j,i) + &
583	kh(k,j,i) * g / prho(k,j,i) * &
584	( prho(k+1,j,i) - prho(k-1,j,i) ) * &
585	dd2zu(k)
586	ENDDO
587	ENDDO
588
589	ELSE
590
591	DO j = nys, nyn
592	DO k = nzb_diff_s_inner(j,i), nzt_diff
593	tend(k,j,i) = tend(k,j,i) - &
594	kh(k,j,i) * g / pt(k,j,i) * &
595	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
596	dd2zu(k)
597	ENDDO
598
599	IF ( use_surface_fluxes ) THEN
600	k = nzb_diff_s_inner(j,i)-1
601	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
602	shf(j,i)
603	ENDIF
604
605	IF ( use_top_fluxes ) THEN
606	k = nzt
607	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
608	tswst(j,i)
609	ENDIF
610	ENDDO
611
612	ENDIF
613
614	ENDIF
615
616	ELSE
617
618	DO j = nys, nyn
619
620	DO k = nzb_diff_s_inner(j,i), nzt_diff
621
622	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
623	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
624	k2 = 0.61_wp * pt(k,j,i)
625	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
626	g / vpt(k,j,i) * &
627	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
628	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
629	) * dd2zu(k)
630	ELSE IF ( cloud_physics ) THEN
631	IF ( ql(k,j,i) == 0.0_wp ) THEN
632	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
633	k2 = 0.61_wp * pt(k,j,i)
634	ELSE
635	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
636	temp = theta * t_d_pt(k)
637	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
638	( q(k,j,i) - ql(k,j,i) ) * &
639	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
640	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
641	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
642	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
643	ENDIF
644	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
645	g / vpt(k,j,i) * &
646	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
647	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
648	) * dd2zu(k)
649	ELSE IF ( cloud_droplets ) THEN
650	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
651	k2 = 0.61_wp * pt(k,j,i)
652	tend(k,j,i) = tend(k,j,i) - &
653	kh(k,j,i) * g / vpt(k,j,i) * &
654	( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &
655	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
656	pt(k,j,i) * ( ql(k+1,j,i) - &
657	ql(k-1,j,i) ) ) * dd2zu(k)
658	ENDIF
659
660	ENDDO
661
662	ENDDO
663
664	IF ( use_surface_fluxes ) THEN
665
666	DO j = nys, nyn
667
668	k = nzb_diff_s_inner(j,i)-1
669
670	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
671	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
672	k2 = 0.61_wp * pt(k,j,i)
673	ELSE IF ( cloud_physics ) THEN
674	IF ( ql(k,j,i) == 0.0_wp ) THEN
675	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
676	k2 = 0.61_wp * pt(k,j,i)
677	ELSE
678	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
679	temp = theta * t_d_pt(k)
680	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
681	( q(k,j,i) - ql(k,j,i) ) * &
682	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
683	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
684	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
685	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
686	ENDIF
687	ELSE IF ( cloud_droplets ) THEN
688	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
689	k2 = 0.61_wp * pt(k,j,i)
690	ENDIF
691
692	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
693	( k1* shf(j,i) + k2 * qsws(j,i) )
694	ENDDO
695
696	ENDIF
697
698	IF ( use_top_fluxes ) THEN
699
700	DO j = nys, nyn
701
702	k = nzt
703
704	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
705	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
706	k2 = 0.61_wp * pt(k,j,i)
707	ELSE IF ( cloud_physics ) THEN
708	IF ( ql(k,j,i) == 0.0_wp ) THEN
709	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
710	k2 = 0.61_wp * pt(k,j,i)
711	ELSE
712	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
713	temp = theta * t_d_pt(k)
714	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
715	( q(k,j,i) - ql(k,j,i) ) * &
716	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
717	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
718	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
719	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
720	ENDIF
721	ELSE IF ( cloud_droplets ) THEN
722	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
723	k2 = 0.61_wp * pt(k,j,i)
724	ENDIF
725
726	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
727	( k1* tswst(j,i) + k2 * qswst(j,i) )
728	ENDDO
729
730	ENDIF
731
732	ENDIF
733
734	ENDIF
735
736	ENDDO
737
738	END SUBROUTINE production_e
739
740
741	!------------------------------------------------------------------------------!
742	! Description:
743	! ------------
744	!> Call for grid point i,j
745	!------------------------------------------------------------------------------!
746	SUBROUTINE production_e_ij( i, j )
747
748	USE arrays_3d, &
749	ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, qsws, qswst, shf, &
750	tend, tswst, u, v, vpt, w
751
752	USE cloud_parameters, &
753	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
754
755	USE control_parameters, &
756	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
757	humidity, kappa, neutral, ocean, pt_reference, &
758	rho_reference, use_single_reference_value, &
759	use_surface_fluxes, use_top_fluxes
760
761	USE grid_variables, &
762	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
763
764	USE indices, &
765	ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, &
766	nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
767
768	IMPLICIT NONE
769
770	INTEGER(iwp) :: i !<
771	INTEGER(iwp) :: j !<
772	INTEGER(iwp) :: k !<
773
774	REAL(wp) :: def !<
775	REAL(wp) :: dudx !<
776	REAL(wp) :: dudy !<
777	REAL(wp) :: dudz !<
778	REAL(wp) :: dvdx !<
779	REAL(wp) :: dvdy !<
780	REAL(wp) :: dvdz !<
781	REAL(wp) :: dwdx !<
782	REAL(wp) :: dwdy !<
783	REAL(wp) :: dwdz !<
784	REAL(wp) :: k1 !<
785	REAL(wp) :: k2 !<
786	REAL(wp) :: km_neutral !<
787	REAL(wp) :: theta !<
788	REAL(wp) :: temp !<
789
790	REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !<
791	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
792	REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !<
793	REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !<
794
795	!
796	!-- Calculate TKE production by shear
797	DO k = nzb_diff_s_outer(j,i), nzt
798
799	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
800	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
801	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
802	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
803	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
804
805	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
806	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
807	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
808	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
809	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
810
811	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
812	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
813	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
814	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
815	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
816
817	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) &
818	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
819	+ 2.0_wp * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
820
821	IF ( def < 0.0_wp ) def = 0.0_wp
822
823	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
824
825	ENDDO
826
827	IF ( constant_flux_layer ) THEN
828
829	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) THEN
830
831	!
832	!-- Position beneath wall
833	!-- (2) - Will allways be executed.
834	!-- 'bottom and wall: use u_0,v_0 and wall functions'
835	k = nzb_diff_s_inner(j,i)-1
836
837	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
838	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
839	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
840	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
841	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
842	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
843	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
844
845	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
846	!
847	!-- Inconsistency removed: as the thermal stratification
848	!-- is not taken into account for the evaluation of the
849	!-- wall fluxes at vertical walls, the eddy viscosity km
850	!-- must not be used for the evaluation of the velocity
851	!-- gradients dudy and dwdy
852	!-- Note: The validity of the new method has not yet
853	!-- been shown, as so far no suitable data for a
854	!-- validation has been available
855	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
856	usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
857	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
858	wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
859	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25_wp &
860	0.5_wp * dy
861	IF ( km_neutral > 0.0_wp ) THEN
862	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
863	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
864	ELSE
865	dudy = 0.0_wp
866	dwdy = 0.0_wp
867	ENDIF
868	ELSE
869	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
870	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
871	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
872	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
873	ENDIF
874
875	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
876	!
877	!-- Inconsistency removed: as the thermal stratification
878	!-- is not taken into account for the evaluation of the
879	!-- wall fluxes at vertical walls, the eddy viscosity km
880	!-- must not be used for the evaluation of the velocity
881	!-- gradients dvdx and dwdx
882	!-- Note: The validity of the new method has not yet
883	!-- been shown, as so far no suitable data for a
884	!-- validation has been available
885	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
886	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
887	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
888	wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
889	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25_wp &
890	0.5_wp * dx
891	IF ( km_neutral > 0.0_wp ) THEN
892	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
893	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
894	ELSE
895	dvdx = 0.0_wp
896	dwdx = 0.0_wp
897	ENDIF
898	ELSE
899	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
900	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
901	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
902	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
903	ENDIF
904
905	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
906	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
907	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
908
909	IF ( def < 0.0_wp ) def = 0.0_wp
910
911	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
912
913	!
914	!-- (3) - will be executed only, if there is at least one level
915	!-- between (2) and (4), i.e. the topography must have a
916	!-- minimum height of 2 dz. Wall fluxes for this case have
917	!-- already been calculated for (2).
918	!-- 'wall only: use wall functions'
919	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
920
921	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
922	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
923	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
924	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
925	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
926	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
927	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
928
929	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
930	!
931	!-- Inconsistency removed: as the thermal stratification
932	!-- is not taken into account for the evaluation of the
933	!-- wall fluxes at vertical walls, the eddy viscosity km
934	!-- must not be used for the evaluation of the velocity
935	!-- gradients dudy and dwdy
936	!-- Note: The validity of the new method has not yet
937	!-- been shown, as so far no suitable data for a
938	!-- validation has been available
939	km_neutral = kappa * ( usvs(k)**2 + &
940	wsvs(k)2 )0.25_wp * 0.5_wp * dy
941	IF ( km_neutral > 0.0_wp ) THEN
942	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
943	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
944	ELSE
945	dudy = 0.0_wp
946	dwdy = 0.0_wp
947	ENDIF
948	ELSE
949	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
950	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
951	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
952	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
953	ENDIF
954
955	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
956	!
957	!-- Inconsistency removed: as the thermal stratification
958	!-- is not taken into account for the evaluation of the
959	!-- wall fluxes at vertical walls, the eddy viscosity km
960	!-- must not be used for the evaluation of the velocity
961	!-- gradients dvdx and dwdx
962	!-- Note: The validity of the new method has not yet
963	!-- been shown, as so far no suitable data for a
964	!-- validation has been available
965	km_neutral = kappa * ( vsus(k)**2 + &
966	wsus(k)2 )0.25_wp * 0.5_wp * dx
967	IF ( km_neutral > 0.0_wp ) THEN
968	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
969	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
970	ELSE
971	dvdx = 0.0_wp
972	dwdx = 0.0_wp
973	ENDIF
974	ELSE
975	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
976	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
977	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
978	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
979	ENDIF
980
981	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
982	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
983	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
984
985	IF ( def < 0.0_wp ) def = 0.0_wp
986
987	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
988
989	ENDDO
990
991	!
992	!-- (4) - will allways be executed.
993	!-- 'special case: free atmosphere' (as for case (0))
994	k = nzb_diff_s_outer(j,i)-1
995
996	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
997	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
998	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
999	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1000	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1001
1002	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1003	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1004	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1005	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1006	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1007
1008	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1009	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1010	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1011	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1012	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1013
1014	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1015	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1016	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1017
1018	IF ( def < 0.0_wp ) def = 0.0_wp
1019
1020	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
1021
1022	ELSE
1023
1024	!
1025	!-- Position without adjacent wall
1026	!-- (1) - will allways be executed.
1027	!-- 'bottom only: use u_0,v_0'
1028	k = nzb_diff_s_inner(j,i)-1
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_0(j,i) - u_0(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_0(j,i) - v_0(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 + dudz2 + dvdz2 &
1050	+ 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	ELSEIF ( use_surface_fluxes ) THEN
1059
1060	k = nzb_diff_s_outer(j,i)-1
1061
1062	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
1063	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
1064	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
1065	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
1066	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
1067
1068	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
1069	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
1070	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
1071	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
1072	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
1073
1074	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
1075	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
1076	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
1077	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
1078	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
1079
1080	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
1081	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
1082	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
1083
1084	IF ( def < 0.0_wp ) def = 0.0_wp
1085
1086	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
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 in
1101	!-- the bottom and top surface layer
1102	DO k = nzb_s_inner(j,i)+1, nzt
1103	tend(k,j,i) = tend(k,j,i) + &
1104	kh(k,j,i) * g / rho_reference * &
1105	( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
1106	ENDDO
1107
1108	ELSE
1109
1110	DO k = nzb_diff_s_inner(j,i), nzt_diff
1111	tend(k,j,i) = tend(k,j,i) - &
1112	kh(k,j,i) * g / pt_reference * &
1113	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
1114	ENDDO
1115
1116	IF ( use_surface_fluxes ) THEN
1117	k = nzb_diff_s_inner(j,i)-1
1118	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
1119	ENDIF
1120
1121	IF ( use_top_fluxes ) THEN
1122	k = nzt
1123	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
1124	ENDIF
1125
1126	ENDIF
1127
1128	ELSE
1129
1130	IF ( ocean ) THEN
1131	!
1132	!-- So far in the ocean no special treatment of density flux in
1133	!-- the bottom and top surface layer
1134	DO k = nzb_s_inner(j,i)+1, nzt
1135	tend(k,j,i) = tend(k,j,i) + &
1136	kh(k,j,i) * g / prho(k,j,i) * &
1137	( prho(k+1,j,i) - prho(k-1,j,i) ) * dd2zu(k)
1138	ENDDO
1139
1140	ELSE
1141
1142	DO k = nzb_diff_s_inner(j,i), nzt_diff
1143	tend(k,j,i) = tend(k,j,i) - &
1144	kh(k,j,i) * g / pt(k,j,i) * &
1145	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
1146	ENDDO
1147
1148	IF ( use_surface_fluxes ) THEN
1149	k = nzb_diff_s_inner(j,i)-1
1150	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
1151	ENDIF
1152
1153	IF ( use_top_fluxes ) THEN
1154	k = nzt
1155	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
1156	ENDIF
1157
1158	ENDIF
1159
1160	ENDIF
1161
1162	ELSE
1163
1164	DO k = nzb_diff_s_inner(j,i), nzt_diff
1165
1166	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1167	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1168	k2 = 0.61_wp * pt(k,j,i)
1169	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
1170	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1171	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
1172	) * dd2zu(k)
1173	ELSE IF ( cloud_physics ) THEN
1174	IF ( ql(k,j,i) == 0.0_wp ) THEN
1175	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1176	k2 = 0.61_wp * pt(k,j,i)
1177	ELSE
1178	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1179	temp = theta * t_d_pt(k)
1180	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1181	( q(k,j,i) - ql(k,j,i) ) * &
1182	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1183	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1184	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1185	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1186	ENDIF
1187	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
1188	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1189	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
1190	) * dd2zu(k)
1191	ELSE IF ( cloud_droplets ) THEN
1192	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1193	k2 = 0.61_wp * pt(k,j,i)
1194	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
1195	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
1196	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
1197	pt(k,j,i) * ( ql(k+1,j,i) - &
1198	ql(k-1,j,i) ) ) * dd2zu(k)
1199	ENDIF
1200	ENDDO
1201
1202	IF ( use_surface_fluxes ) THEN
1203	k = nzb_diff_s_inner(j,i)-1
1204
1205	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1206	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1207	k2 = 0.61_wp * pt(k,j,i)
1208	ELSE IF ( cloud_physics ) THEN
1209	IF ( ql(k,j,i) == 0.0_wp ) THEN
1210	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1211	k2 = 0.61_wp * pt(k,j,i)
1212	ELSE
1213	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1214	temp = theta * t_d_pt(k)
1215	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1216	( q(k,j,i) - ql(k,j,i) ) * &
1217	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1218	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1219	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1220	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1221	ENDIF
1222	ELSE IF ( cloud_droplets ) THEN
1223	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1224	k2 = 0.61_wp * pt(k,j,i)
1225	ENDIF
1226
1227	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1228	( k1* shf(j,i) + k2 * qsws(j,i) )
1229	ENDIF
1230
1231	IF ( use_top_fluxes ) THEN
1232	k = nzt
1233
1234	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
1235	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1236	k2 = 0.61_wp * pt(k,j,i)
1237	ELSE IF ( cloud_physics ) THEN
1238	IF ( ql(k,j,i) == 0.0_wp ) THEN
1239	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
1240	k2 = 0.61_wp * pt(k,j,i)
1241	ELSE
1242	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
1243	temp = theta * t_d_pt(k)
1244	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
1245	( q(k,j,i) - ql(k,j,i) ) * &
1246	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
1247	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
1248	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
1249	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
1250	ENDIF
1251	ELSE IF ( cloud_droplets ) THEN
1252	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
1253	k2 = 0.61_wp * pt(k,j,i)
1254	ENDIF
1255
1256	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
1257	( k1* tswst(j,i) + k2 * qswst(j,i) )
1258	ENDIF
1259
1260	ENDIF
1261
1262	ENDIF
1263
1264	END SUBROUTINE production_e_ij
1265
1266
1267	!------------------------------------------------------------------------------!
1268	! Description:
1269	! ------------
1270	!> @todo Missing subroutine description.
1271	!------------------------------------------------------------------------------!
1272	SUBROUTINE production_e_init
1273
1274	USE arrays_3d, &
1275	ONLY: kh, km, u, us, usws, v, vsws, zu
1276
1277	USE control_parameters, &
1278	ONLY: constant_flux_layer, kappa
1279
1280	USE indices, &
1281	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner, &
1282	nzb_v_inner
1283
1284	IMPLICIT NONE
1285
1286	INTEGER(iwp) :: i !<
1287	INTEGER(iwp) :: j !<
1288	INTEGER(iwp) :: ku !<
1289	INTEGER(iwp) :: kv !<
1290
1291	IF ( constant_flux_layer ) THEN
1292
1293	IF ( first_call ) THEN
1294	ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
1295	u_0 = 0.0_wp ! just to avoid access of uninitialized memory
1296	v_0 = 0.0_wp ! within exchange_horiz_2d
1297	first_call = .FALSE.
1298	ENDIF
1299
1300	!
1301	!-- Calculate a virtual velocity at the surface in a way that the
1302	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
1303	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
1304	!-- production term at k=1 (see production_e_ij).
1305	!-- The velocity gradient has to be limited in case of too small km
1306	!-- (otherwise the timestep may be significantly reduced by large
1307	!-- surface winds).
1308	!-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are
1309	!-- not available in case of non-cyclic boundary conditions.
1310	!-- WARNING: the exact analytical solution would require the determination
1311	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
1312	!$OMP PARALLEL DO PRIVATE( ku, kv )
1313	DO i = nxl, nxr+1
1314	DO j = nys, nyn+1
1315
1316	ku = nzb_u_inner(j,i)+1
1317	kv = nzb_v_inner(j,i)+1
1318
1319	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
1320	( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) + &
1321	1.0E-20_wp )
1322	! ( us(j,i) * kappa * zu(1) )
1323	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
1324	( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) + &
1325	1.0E-20_wp )
1326	! ( us(j,i) * kappa * zu(1) )
1327
1328	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
1329	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
1330	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
1331	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
1332
1333	ENDDO
1334	ENDDO
1335
1336	CALL exchange_horiz_2d( u_0 )
1337	CALL exchange_horiz_2d( v_0 )
1338
1339	ENDIF
1340
1341	END SUBROUTINE production_e_init
1342
1343	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |