Home

Context Navigation

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

Last change on this file since 2100 was 2032, checked in by knoop, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 89.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |