Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |