Home

Context Navigation

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

Last change on this file since 1370 was 1354, checked in by heinze, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 88.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |