Home

Context Navigation

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

Last change on this file since 1949 was 1874, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 89.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |