Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |