Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |