Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |