Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |