Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |