Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |