Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |