Home

Context Navigation

source: palm/trunk/SOURCE/production_e_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 89.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |