Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |