Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |