Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |