Home

Context Navigation

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

Last change on this file since 2335 was 2329, checked in by knoop, 8 years ago
Bugfix for topography usage with anelastic approximation and boussinesq approximation with air density != 1
Property svn:keywords set to `Id`
File size: 85.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |