Home

Context Navigation

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

Last change on this file since 1350 was 1343, checked in by kanani, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 88.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |