Home

Context Navigation

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

Last change on this file since 2282 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 82.6 KB

Rev	Line
[1873]	1	!> @file production_e.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	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.
[1036]	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	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[2233]	22	!
	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: production_e.f90 2233 2017-05-30 18:08:54Z schwenkel $
	27	!
[2233]	28	! 2232 2017-05-30 17:47:52Z suehring
	29	! Adjustments to new surface concept
	30	!
[2127]	31	! 2126 2017-01-20 15:54:21Z raasch
	32	! density in ocean case replaced by potential density
	33	!
[2119]	34	! 2118 2017-01-17 16:38:49Z raasch
	35	! OpenACC version of subroutine removed
	36	!
[2032]	37	! 2031 2016-10-21 15:11:58Z knoop
	38	! renamed variable rho to rho_ocean
	39	!
[2001]	40	! 2000 2016-08-20 18:09:15Z knoop
	41	! Forced header and separation lines into 80 columns
	42	!
[1874]	43	! 1873 2016-04-18 14:50:06Z maronga
	44	! Module renamed (removed _mod)
	45	!
	46	!
[1851]	47	! 1850 2016-04-08 13:29:27Z maronga
	48	! Module renamed
	49	!
	50	!
[1692]	51	! 1691 2015-10-26 16:17:44Z maronga
	52	! Renamed prandtl_layer to constant_flux_layer.
	53	!
[1683]	54	! 1682 2015-10-07 23:56:08Z knoop
	55	! Code annotations made doxygen readable
	56	!
[1375]	57	! 1374 2014-04-25 12:55:07Z raasch
	58	! nzb_s_outer removed from acc-present-list
	59	!
[1354]	60	! 1353 2014-04-08 15:21:23Z heinze
	61	! REAL constants provided with KIND-attribute
	62	!
[1343]	63	! 1342 2014-03-26 17:04:47Z kanani
	64	! REAL constants defined as wp-kind
	65	!
[1321]	66	! 1320 2014-03-20 08:40:49Z raasch
[1320]	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
[110]	74	!
[1258]	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	!
[1182]	79	! 1179 2013-06-14 05:57:58Z raasch
	80	! use_reference renamed use_single_reference_value
	81	!
[1132]	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	!
[1037]	86	! 1036 2012-10-22 13:43:42Z raasch
	87	! code put under GPL (PALM 3.9)
	88	!
[1017]	89	! 1015 2012-09-27 09:23:24Z raasch
	90	! accelerator version (*_acc) added
	91	!
[1008]	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	!
[941]	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	!
[1]	100	! Revision 1.1 1997/09/19 07:45:35 raasch
	101	! Initial revision
	102	!
	103	!
	104	! Description:
	105	! ------------
[1682]	106	!> Production terms (shear + buoyancy) of the TKE.
[1691]	107	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
[1682]	108	!> not considered well!
[1]	109	!------------------------------------------------------------------------------!
[1682]	110	MODULE production_e_mod
[1]	111
[1320]	112	USE kinds
	113
[1]	114	PRIVATE
[2118]	115	PUBLIC production_e, production_e_init
[56]	116
[1]	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	!------------------------------------------------------------------------------!
[1682]	130	! Description:
	131	! ------------
	132	!> Call for all grid points
[1]	133	!------------------------------------------------------------------------------!
	134	SUBROUTINE production_e
	135
[1320]	136	USE arrays_3d, &
[2232]	137	ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
[1]	138
[1320]	139	USE cloud_parameters, &
	140	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
	141
	142	USE control_parameters, &
[1691]	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
[1320]	147
	148	USE grid_variables, &
[2232]	149	ONLY: ddx, dx, ddy, dy
[1320]	150
	151	USE indices, &
[2232]	152	ONLY: nxl, nxr, nys, nyn, nzb, nzt, wall_flags_0
[1320]	153
[2232]	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
[1]	158	IMPLICIT NONE
	159
[2232]	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
[1]	167
[2232]	168	REAL(wp) :: def !<
	169	REAL(wp) :: flag !< flag to mask topography
[1682]	170	REAL(wp) :: k1 !<
	171	REAL(wp) :: k2 !<
[2232]	172	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
[1682]	173	REAL(wp) :: theta !<
	174	REAL(wp) :: temp !<
[2232]	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"
[1]	180
[2232]	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
[1]	190
[2232]	191	DO i = nxl, nxr
[53]	192
[2232]	193	IF ( constant_flux_layer ) THEN
[940]	194
	195	!
[2232]	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
[1]	201
[2232]	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)
[1]	215
[2232]	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)
[1]	221
[2232]	222	ENDDO
[1]	223	ENDDO
	224
	225	!
[55]	226	!-- Position beneath wall
	227	!-- (2) - Will allways be executed.
	228	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	229	DO j = nys, nyn
[1007]	230	!
[2232]	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
[1007]	244	!
[2232]	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
[1]	251	!
[2232]	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)
[1]	297
[2232]	298	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	299	* 0.5_wp * dx
[1007]	300	!
[2232]	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)
[1]	315
[2232]	316	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	317	* 0.5_wp * dx
[1007]	318	!
[2232]	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)
[1]	333
[2232]	334	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	335	* 0.5_wp * dx
[1]	336	!
[2232]	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)
[1]	407
[2232]	408	ENDDO
[1]	409
	410	ENDDO
	411
	412	DO j = nys, nyn
[2232]	413	DO k = nzb+1, nzt
[1]	414
[2232]	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) )
[1]	420
[1342]	421	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	422
[2232]	423	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1007]	424
[2232]	425	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
[1]	426
[2232]	427	ENDDO
[1]	428	ENDDO
	429
[37]	430	ELSEIF ( use_surface_fluxes ) THEN
	431
	432	DO j = nys, nyn
[2232]	433	!
	434	!-- Calculate TKE production by shear. Here, no additional
	435	!-- wall-bounded code is considered.
	436	!-- Why?
	437	DO k = nzb+1, nzt
[37]	438
[2232]	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)
[37]	445
[2232]	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)
[37]	452
[2232]	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	)
[37]	469
[2232]	470	IF ( def < 0.0_wp ) def = 0.0_wp
[37]	471
[2232]	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
[37]	477	ENDDO
	478
[1]	479	ENDIF
	480
	481	!
[940]	482	!-- If required, calculate TKE production by buoyancy
	483	IF ( .NOT. neutral ) THEN
[1]	484
[940]	485	IF ( .NOT. humidity ) THEN
[1]	486
[1179]	487	IF ( use_single_reference_value ) THEN
[940]	488
	489	IF ( ocean ) THEN
[97]	490	!
[940]	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
[2232]	494	DO k = nzb+1, nzt
[2126]	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) ) * &
[2232]	498	dd2zu(k) * &
	499	MERGE( 1.0_wp, 0.0_wp, &
	500	BTEST( wall_flags_0(k,j,i), 0 ) &
	501	)
[940]	502	ENDDO
[97]	503	ENDDO
	504
[940]	505	ELSE
[97]	506
[940]	507	DO j = nys, nyn
[2232]	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	)
[940]	522	ENDDO
[97]	523
[940]	524	IF ( use_surface_fluxes ) THEN
[2232]	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
[940]	554	ENDIF
[97]	555
[940]	556	IF ( use_top_fluxes ) THEN
[2232]	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
[940]	564	ENDIF
	565	ENDDO
[57]	566
[940]	567	ENDIF
[57]	568
[940]	569	ELSE
[1]	570
[940]	571	IF ( ocean ) THEN
[97]	572	!
[940]	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
[2232]	576	DO k = nzb+1, nzt
[2126]	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) ) * &
[2232]	580	dd2zu(k) * &
	581	MERGE( 1.0_wp, 0.0_wp, &
	582	BTEST( wall_flags_0(k,j,i), 0 ) &
	583	)
[940]	584	ENDDO
[97]	585	ENDDO
	586
[940]	587	ELSE
[97]	588
[940]	589	DO j = nys, nyn
[2232]	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	)
[940]	604	ENDDO
	605
	606	IF ( use_surface_fluxes ) THEN
[2232]	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
[940]	636	ENDIF
	637
	638	IF ( use_top_fluxes ) THEN
[2232]	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
[940]	646	ENDIF
[97]	647	ENDDO
	648
[940]	649	ENDIF
[97]	650
	651	ENDIF
[1]	652
[940]	653	ELSE
[57]	654
[940]	655	DO j = nys, nyn
[1]	656
[2232]	657	DO k = nzb+1, nzt
	658	!
	659	!-- Flag 9 is used to mask top fluxes, flag 30 to mask
	660	!-- surface fluxes
[1007]	661	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	662	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	663	k2 = 0.61_wp * pt(k,j,i)
[1007]	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) ) &
[2232]	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	)
[1007]	675	ELSE IF ( cloud_physics ) THEN
[1342]	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)
[940]	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)
[1342]	682	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
[2232]	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 * &
[940]	686	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	687	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	688	ENDIF
[1007]	689	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	690	g / vpt(k,j,i) * &
[940]	691	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	692	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
[2232]	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	)
[1007]	700	ELSE IF ( cloud_droplets ) THEN
[1342]	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)
[2232]	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	)
[1007]	715	ENDIF
	716
[940]	717	ENDDO
	718
[1]	719	ENDDO
	720
[940]	721	IF ( use_surface_fluxes ) THEN
[1]	722
[940]	723	DO j = nys, nyn
[2232]	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)
[1]	731
[2232]	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
[1]	753
[2232]	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
[1342]	767	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	768	k2 = 0.61_wp * pt(k,j,i)
[2232]	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)
[940]	786	ENDIF
[1]	787
[2232]	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
[940]	826	ENDDO
[1]	827
[940]	828	ENDIF
[1]	829
[940]	830	IF ( use_top_fluxes ) THEN
[19]	831
[940]	832	DO j = nys, nyn
[19]	833
[2232]	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)
[19]	838
[2232]	839	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	840	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	841	k2 = 0.61_wp * pt(k,j,i)
[2232]	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 * &
[1353]	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 * &
[940]	853	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[2232]	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)
[940]	859	ENDIF
[19]	860
[2232]	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
[940]	867	ENDDO
[19]	868
[940]	869	ENDIF
	870
[19]	871	ENDIF
	872
[1]	873	ENDIF
	874
	875	ENDDO
	876
	877	END SUBROUTINE production_e
	878
	879
	880	!------------------------------------------------------------------------------!
[1682]	881	! Description:
	882	! ------------
	883	!> Call for grid point i,j
[1]	884	!------------------------------------------------------------------------------!
	885	SUBROUTINE production_e_ij( i, j )
	886
[1320]	887	USE arrays_3d, &
[2232]	888	ONLY: ddzw, dd2zu, kh, km, prho, pt, q, ql, tend, u, v, vpt, w
[449]	889
[1320]	890	USE cloud_parameters, &
	891	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
	892
	893	USE control_parameters, &
[1691]	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
[1320]	898
	899	USE grid_variables, &
[2232]	900	ONLY: ddx, dx, ddy, dy
[1320]	901
	902	USE indices, &
[2232]	903	ONLY: nxl, nxr, nys, nyn, nzb, nzb, nzt, wall_flags_0
[1320]	904
[2232]	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
[1]	909	IMPLICIT NONE
	910
[2232]	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
[1]	918
[1682]	919	REAL(wp) :: def !<
[2232]	920	REAL(wp) :: flag !< flag to mask topography
[1682]	921	REAL(wp) :: k1 !<
	922	REAL(wp) :: k2 !<
[2232]	923	REAL(wp) :: km_neutral !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
[1682]	924	REAL(wp) :: theta !<
	925	REAL(wp) :: temp !<
[2232]	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"
[1]	931
[53]	932
[2232]	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
[1]	942
	943
[2232]	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
[1]	950
[2232]	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)
[1]	956
[2232]	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)
[1]	962
[2232]	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)
[1]	968
[2232]	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)
[1007]	982
[2232]	983	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	984	* 0.5_wp * dy
[1]	985	!
[2232]	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)
[1]	1000
[2232]	1001	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	1002	* 0.5_wp * dy
[1007]	1003	!
[2232]	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)
[1]	1018
[2232]	1019	km_neutral = kappa * ( usvs2 + wsvs2 )**0.25_wp &
	1020	* 0.5_wp * dy
[1007]	1021	!
[2232]	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)
[1]	1040
[2232]	1041	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	1042	* 0.5_wp * dx
[1]	1043	!
[2232]	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)
[1]	1058
[2232]	1059	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	1060	* 0.5_wp * dx
[1007]	1061	!
[2232]	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)
[1]	1076
[2232]	1077	km_neutral = kappa * ( vsus2 + wsus2 )**0.25_wp &
	1078	* 0.5_wp * dx
[1007]	1079	!
[2232]	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)
[1]	1102
[2232]	1103	ENDDO
[1]	1104	!
[2232]	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)
[1]	1149
[2232]	1150	ENDDO
[1]	1151
[2232]	1152	DO k = nzb+1, nzt
[1]	1153
[2232]	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) )
[1]	1158
[1342]	1159	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	1160
[2232]	1161	flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1]	1162
[2232]	1163	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag
[1]	1164
[2232]	1165	ENDDO
	1166
	1167	ELSEIF ( use_surface_fluxes ) THEN
[1]	1168	!
[2232]	1169	!-- Calculate TKE production by shear. Here, no additional
	1170	!-- wall-bounded code is considered.
	1171	!-- Why?
	1172	DO k = nzb+1, nzt
[1]	1173
[2232]	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)
[1]	1179
[2232]	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)
[1]	1185
[2232]	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)
[1]	1191
[2232]	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) )
[1]	1196
[1342]	1197	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	1198
[2232]	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
[1]	1202
[2232]	1203	ENDDO
[1]	1204
	1205	ENDIF
	1206
	1207	!
[940]	1208	!-- If required, calculate TKE production by buoyancy
	1209	IF ( .NOT. neutral ) THEN
[1]	1210
[940]	1211	IF ( .NOT. humidity ) THEN
[19]	1212
[1179]	1213	IF ( use_single_reference_value ) THEN
[940]	1214
	1215	IF ( ocean ) THEN
[97]	1216	!
[940]	1217	!-- So far in the ocean no special treatment of density flux in
	1218	!-- the bottom and top surface layer
[2232]	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	)
[940]	1227	ENDDO
[97]	1228
[940]	1229	ELSE
[97]	1230
[2232]	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
[940]	1245	ENDDO
[1]	1246
[940]	1247	IF ( use_surface_fluxes ) THEN
[2232]	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
[940]	1277	ENDIF
[19]	1278
[940]	1279	IF ( use_top_fluxes ) THEN
[2232]	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
[940]	1287	ENDIF
	1288
[97]	1289	ENDIF
	1290
[940]	1291	ELSE
[57]	1292
[940]	1293	IF ( ocean ) THEN
[97]	1294	!
[940]	1295	!-- So far in the ocean no special treatment of density flux in
	1296	!-- the bottom and top surface layer
[2232]	1297	DO k = nzb+1, nzt
[2126]	1298	tend(k,j,i) = tend(k,j,i) + &
[2232]	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	)
[940]	1305	ENDDO
[97]	1306
[940]	1307	ELSE
[97]	1308
[2232]	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	)
[940]	1322	ENDDO
[57]	1323
[940]	1324	IF ( use_surface_fluxes ) THEN
[2232]	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
[940]	1354	ENDIF
[57]	1355
[940]	1356	IF ( use_top_fluxes ) THEN
[2232]	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
[940]	1364	ENDIF
	1365
[97]	1366	ENDIF
	1367
[57]	1368	ENDIF
	1369
[940]	1370	ELSE
[57]	1371
[2232]	1372	DO k = nzb+1, nzt
	1373	!
	1374	!-- Flag 9 is used to mask top fluxes, flag 30 to mask surface fluxes
[1007]	1375	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1376	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1377	k2 = 0.61_wp * pt(k,j,i)
[2232]	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
[1007]	1389	ELSE IF ( cloud_physics ) THEN
[1342]	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)
[940]	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)
[2232]	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 * &
[940]	1400	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	1401	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	1402	ENDIF
[2232]	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	)
[1007]	1413	ELSE IF ( cloud_droplets ) THEN
[1342]	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)
[2232]	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	)
[1007]	1427	ENDIF
[940]	1428	ENDDO
[19]	1429
[940]	1430	IF ( use_surface_fluxes ) THEN
[2232]	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)
[1]	1438
[2232]	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)
[1342]	1492	k2 = 0.61_wp * pt(k,j,i)
[940]	1493	ENDIF
	1494
[2232]	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
[1]	1533	ENDIF
	1534
[940]	1535	IF ( use_top_fluxes ) THEN
[2232]	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)
[1]	1540
[2232]	1541
	1542
	1543	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1544	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1545	k2 = 0.61_wp * pt(k,j,i)
[2232]	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 * &
[940]	1557	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[2232]	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)
[940]	1563	ENDIF
	1564
[2232]	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
[19]	1570	ENDIF
	1571
	1572	ENDIF
	1573
[1]	1574	ENDIF
	1575
	1576	END SUBROUTINE production_e_ij
	1577
	1578
[1682]	1579	!------------------------------------------------------------------------------!
	1580	! Description:
	1581	! ------------
	1582	!> @todo Missing subroutine description.
	1583	!------------------------------------------------------------------------------!
[1]	1584	SUBROUTINE production_e_init
	1585
[1320]	1586	USE arrays_3d, &
[2232]	1587	ONLY: kh, km, u, v, zu
[1]	1588
[1320]	1589	USE control_parameters, &
[1691]	1590	ONLY: constant_flux_layer, kappa
[1320]	1591
	1592	USE indices, &
	1593	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner, &
	1594	nzb_v_inner
	1595
[2232]	1596	USE surface_mod, &
	1597	ONLY : surf_def_h, surf_def_v, surf_lsm_h, surf_usm_h
	1598
[1]	1599	IMPLICIT NONE
	1600
[2232]	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
[1]	1606
[1691]	1607	IF ( constant_flux_layer ) THEN
[1]	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).
[106]	1616	!-- not available in case of non-cyclic boundary conditions.
[1]	1617	!-- WARNING: the exact analytical solution would require the determination
	1618	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
[2232]	1619	!-- Default surfaces, upward-facing
	1620	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
	1621	DO m = 1, surf_def_h(0)%ns
[1]	1622
[2232]	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 )
[1]	1637
[2232]	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)
[1]	1641
[2232]	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)
[1]	1645
	1646	ENDDO
[2232]	1647	!
	1648	!-- Default surfaces, downward-facing
	1649	!$OMP PARALLEL DO PRIVATE(i,j,k,m)
	1650	DO m = 1, surf_def_h(1)%ns
[1]	1651
[2232]	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 )
[1]	1668
[2232]	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
[1]	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 |