Home

Context Navigation

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

Last change on this file since 1203 was 1182, checked in by raasch, 11 years ago
last commit documented, rc-file for example run updated
Property svn:keywords set to `Id`
File size: 83.2 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	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[110]	22	!
[1182]	23	!
[110]	24	! Former revisions:
	25	! -----------------
	26	! $Id: production_e.f90 1182 2013-06-14 09:07:24Z witha $
	27	!
[1182]	28	! 1179 2013-06-14 05:57:58Z raasch
	29	! use_reference renamed use_single_reference_value
	30	!
[1132]	31	! 1128 2013-04-12 06:19:32Z raasch
	32	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	33	! j_north
	34	!
[1037]	35	! 1036 2012-10-22 13:43:42Z raasch
	36	! code put under GPL (PALM 3.9)
	37	!
[1017]	38	! 1015 2012-09-27 09:23:24Z raasch
	39	! accelerator version (*_acc) added
	40	!
[1008]	41	! 1007 2012-09-19 14:30:36Z franke
	42	! Bugfix: calculation of buoyancy production has to consider the liquid water
	43	! mixing ratio in case of cloud droplets
	44	!
[941]	45	! 940 2012-07-09 14:31:00Z raasch
	46	! TKE production by buoyancy can be switched off in case of runs with pure
	47	! neutral stratification
	48	!
[760]	49	! 759 2011-09-15 13:58:31Z raasch
	50	! initialization of u_0, v_0
	51	!
[668]	52	! 667 2010-12-23 12:06:00Z suehring/gryschka
	53	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
	54	!
[482]	55	! 449 2010-02-02 11:23:59Z raasch
	56	! test output from rev 410 removed
	57	!
[392]	58	! 388 2009-09-23 09:40:33Z raasch
	59	! Bugfix: wrong sign in buoyancy production of ocean part in case of not using
	60	! the reference density (only in 3D routine production_e)
	61	! Bugfix to avoid zero division by km_neutral
	62	!
[226]	63	! 208 2008-10-20 06:02:59Z raasch
	64	! Bugfix concerning the calculation of velocity gradients at vertical walls
	65	! in case of diabatic conditions
	66	!
[198]	67	! 187 2008-08-06 16:25:09Z letzel
	68	! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
	69	! consistency with subroutine wall_fluxes
	70	!
[139]	71	! 124 2007-10-19 15:47:46Z raasch
	72	! Bugfix: calculation of density flux in the ocean now starts from nzb+1
	73	!
[110]	74	! 108 2007-08-24 15:10:38Z letzel
[106]	75	! Bugfix: wrong sign removed from the buoyancy production term in the case
	76	! use_reference = .T.,
	77	! u_0 and v_0 are calculated for nxr+1, nyn+1 also (otherwise these values are
	78	! not available in case of non-cyclic boundary conditions)
[108]	79	! Bugfix for ocean density flux at bottom
[39]	80	!
[98]	81	! 97 2007-06-21 08:23:15Z raasch
	82	! Energy production by density flux (in ocean) added
	83	! use_pt_reference renamed use_reference
	84	!
[77]	85	! 75 2007-03-22 09:54:05Z raasch
	86	! Wall functions now include diabatic conditions, call of routine wall_fluxes_e,
	87	! reference temperature pt_reference can be used in buoyancy term,
	88	! moisture renamed humidity
	89	!
[39]	90	! 37 2007-03-01 08:33:54Z raasch
[19]	91	! Calculation extended for gridpoint nzt, extended for given temperature /
[37]	92	! humidity fluxes at the top, wall-part is now executed in case that a
	93	! Prandtl-layer is switched on (instead of surfaces fluxes switched on)
[1]	94	!
[3]	95	! RCS Log replace by Id keyword, revision history cleaned up
	96	!
[1]	97	! Revision 1.21 2006/04/26 12:45:35 raasch
	98	! OpenMP parallelization of production_e_init
	99	!
	100	! Revision 1.1 1997/09/19 07:45:35 raasch
	101	! Initial revision
	102	!
	103	!
	104	! Description:
	105	! ------------
	106	! Production terms (shear + buoyancy) of the TKE
[37]	107	! WARNING: the case with prandtl_layer = F and use_surface_fluxes = T is
	108	! not considered well!
[1]	109	!------------------------------------------------------------------------------!
	110
[56]	111	USE wall_fluxes_mod
	112
[1]	113	PRIVATE
[1015]	114	PUBLIC production_e, production_e_acc, production_e_init
[56]	115
[1]	116	LOGICAL, SAVE :: first_call = .TRUE.
	117
	118	REAL, DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0, v_0
	119
	120	INTERFACE production_e
	121	MODULE PROCEDURE production_e
	122	MODULE PROCEDURE production_e_ij
	123	END INTERFACE production_e
	124
[1015]	125	INTERFACE production_e_acc
	126	MODULE PROCEDURE production_e_acc
	127	END INTERFACE production_e_acc
	128
[1]	129	INTERFACE production_e_init
	130	MODULE PROCEDURE production_e_init
	131	END INTERFACE production_e_init
	132
	133	CONTAINS
	134
	135
	136	!------------------------------------------------------------------------------!
	137	! Call for all grid points
	138	!------------------------------------------------------------------------------!
	139	SUBROUTINE production_e
	140
	141	USE arrays_3d
	142	USE cloud_parameters
	143	USE control_parameters
	144	USE grid_variables
	145	USE indices
	146	USE statistics
	147
	148	IMPLICIT NONE
	149
	150	INTEGER :: i, j, k
	151
	152	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
[208]	153	k1, k2, km_neutral, theta, temp
[1]	154
[56]	155	! REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs
	156	REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, 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
	164	! CALL wall_fluxes_e( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
	165	! CALL wall_fluxes_e( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
	166	! CALL wall_fluxes_e( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
	167	! CALL wall_fluxes_e( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
	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
	178	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	179	dudy = 0.25 * ( 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 * ( 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)
	183
	184	dvdx = 0.25 * ( 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 * ( 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)
	189
	190	dwdx = 0.25 * ( 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 * ( 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)
	195
	196	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	197	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	198	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	199
	200	IF ( def < 0.0 ) def = 0.0
	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
	215	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	216	THEN
	217
	218	k = nzb_diff_s_inner(j,i) - 1
	219	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	220	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	221	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	222	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	223	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	224	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	225	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	226
[1]	227	IF ( wall_e_y(j,i) /= 0.0 ) 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, &
	238	usvs, 1.0, 0.0, 0.0, 0.0 )
	239	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	240	wsvs, 0.0, 0.0, 1.0, 0.0 )
[208]	241	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25 &
	242	0.5 * dy
[364]	243	IF ( km_neutral > 0.0 ) THEN
	244	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	245	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	246	ELSE
	247	dudy = 0.0
	248	dwdy = 0.0
	249	ENDIF
[1]	250	ELSE
	251	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	252	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	253	dwdy = 0.25 * ( 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
	257	IF ( wall_e_x(j,i) /= 0.0 ) 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, &
	268	vsus, 0.0, 1.0, 0.0, 0.0 )
	269	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	270	wsus, 0.0, 0.0, 0.0, 1.0 )
[208]	271	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25 &
	272	0.5 * dx
[364]	273	IF ( km_neutral > 0.0 ) THEN
	274	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	275	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	276	ELSE
	277	dvdx = 0.0
	278	dwdx = 0.0
	279	ENDIF
[1]	280	ELSE
	281	dvdx = 0.25 * ( 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 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	284	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	285	ENDIF
	286
	287	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	288	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	289	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	290
	291	IF ( def < 0.0 ) def = 0.0
	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
[53]	306	dudz = 0.5 * ( 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 * ( 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)
	311	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	312
[1]	313	IF ( wall_e_y(j,i) /= 0.0 ) 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 + &
	324	wsvs(k)2 )0.25 * 0.5 * dy
[364]	325	IF ( km_neutral > 0.0 ) THEN
	326	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	327	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	328	ELSE
	329	dudy = 0.0
	330	dwdy = 0.0
	331	ENDIF
[1]	332	ELSE
	333	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	334	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	335	dwdy = 0.25 * ( 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
	339	IF ( wall_e_x(j,i) /= 0.0 ) 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 + &
	350	wsus(k)2 )0.25 * 0.5 * dx
[364]	351	IF ( km_neutral > 0.0 ) THEN
	352	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	353	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	354	ELSE
	355	dvdx = 0.0
	356	dwdx = 0.0
	357	ENDIF
[1]	358	ELSE
	359	dvdx = 0.25 * ( 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 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	362	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	363	ENDIF
	364
	365	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	366	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	367	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	368
	369	IF ( def < 0.0 ) def = 0.0
	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
	384	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) &
	385	THEN
	386
	387	k = nzb_diff_s_outer(j,i)-1
	388
	389	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	390	dudy = 0.25 * ( 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 * ( 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)
	394
	395	dvdx = 0.25 * ( 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 * ( 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)
	400
	401	dwdx = 0.25 * ( 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 * ( 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)
	406
	407	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	408	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	409	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	410
	411	IF ( def < 0.0 ) def = 0.0
	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
	425	IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
	426	THEN
	427
	428	k = nzb_diff_s_inner(j,i)-1
	429
	430	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	431	dudy = 0.25 * ( 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 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	434	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	435
	436	dvdx = 0.25 * ( 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 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	440	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	441
	442	dwdx = 0.25 * ( 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 * ( 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)
	447
	448	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	449	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	450	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	451
	452	IF ( def < 0.0 ) def = 0.0
	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
	466	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	467	dudy = 0.25 * ( 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 * ( 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)
	471
	472	dvdx = 0.25 * ( 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 * ( 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)
	477
	478	dwdx = 0.25 * ( 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 * ( 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)
	483
	484	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	485	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	486	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	487
	488	IF ( def < 0.0 ) def = 0.0
	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
[1]	591	k1 = 1.0 + 0.61 * q(k,j,i)
	592	k2 = 0.61 * 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
[940]	599	IF ( ql(k,j,i) == 0.0 ) THEN
	600	k1 = 1.0 + 0.61 * q(k,j,i)
	601	k2 = 0.61 * pt(k,j,i)
	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)
	605	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	606	( q(k,j,i) - ql(k,j,i) ) * &
	607	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	608	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	609	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	610	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	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
	618	k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	619	k2 = 0.61 * pt(k,j,i)
	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
[1]	639	k1 = 1.0 + 0.61 * q(k,j,i)
	640	k2 = 0.61 * pt(k,j,i)
[1007]	641	ELSE IF ( cloud_physics ) THEN
[940]	642	IF ( ql(k,j,i) == 0.0 ) THEN
	643	k1 = 1.0 + 0.61 * q(k,j,i)
	644	k2 = 0.61 * pt(k,j,i)
	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)
	648	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	649	( q(k,j,i) - ql(k,j,i) ) * &
	650	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	651	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	652	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	653	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	654	ENDIF
[1007]	655	ELSE IF ( cloud_droplets ) THEN
	656	k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	657	k2 = 0.61 * 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
[19]	673	k1 = 1.0 + 0.61 * q(k,j,i)
	674	k2 = 0.61 * pt(k,j,i)
[1007]	675	ELSE IF ( cloud_physics ) THEN
[940]	676	IF ( ql(k,j,i) == 0.0 ) THEN
	677	k1 = 1.0 + 0.61 * q(k,j,i)
	678	k2 = 0.61 * pt(k,j,i)
	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)
	682	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	683	( q(k,j,i) - ql(k,j,i) ) * &
	684	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	685	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	686	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	687	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	688	ENDIF
[1007]	689	ELSE IF ( cloud_droplets ) THEN
	690	k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	691	k2 = 0.61 * 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
	714	USE arrays_3d
	715	USE cloud_parameters
	716	USE control_parameters
	717	USE grid_variables
	718	USE indices
	719	USE statistics
	720
	721	IMPLICIT NONE
	722
	723	INTEGER :: i, j, k
	724
	725	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
	726	k1, k2, km_neutral, theta, temp
	727
	728	!$acc declare create ( usvs, vsus, wsus, wsvs )
	729	REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs
	730	! REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs
	731
	732	!
	733	!-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
	734	!-- vertical walls, if neccessary
	735	!-- CAUTION: results are slightly different from the ij-version!!
	736	!-- ij-version should be called further below within the ij-loops!!
	737	IF ( topography /= 'flat' ) THEN
	738	CALL wall_fluxes_e_acc( usvs, 1.0, 0.0, 0.0, 0.0, wall_e_y )
	739	CALL wall_fluxes_e_acc( wsvs, 0.0, 0.0, 1.0, 0.0, wall_e_y )
	740	CALL wall_fluxes_e_acc( vsus, 0.0, 1.0, 0.0, 0.0, wall_e_x )
	741	CALL wall_fluxes_e_acc( wsus, 0.0, 0.0, 0.0, 1.0, wall_e_x )
	742	ENDIF
	743
	744
	745	!
	746	!-- Calculate TKE production by shear
	747	!$acc kernels present( ddzw, dd2zu, kh, km, nzb_diff_s_inner, nzb_diff_s_outer ) &
	748	!$acc present( nzb_s_inner, nzb_s_outer, pt, q, ql, qsws, qswst, rho ) &
	749	!$acc present( shf, tend, tswst, u, v, vpt, w, wall_e_x, wall_e_y ) &
	750	!$acc copyin( u_0, v_0 )
	751	!$acc loop
[1128]	752	DO i = i_left, i_right
	753	DO j = j_south, j_north
[1015]	754	!$acc loop vector( 32 )
	755	DO k = 1, nzt
	756
	757	IF ( k >= nzb_diff_s_outer(j,i) ) THEN
	758
	759	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	760	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	761	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	762	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	763	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	764
	765	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	766	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	767	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	768	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	769	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	770
	771	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	772	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	773	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	774	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	775	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	776
	777	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	778	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	779	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	780
	781	IF ( def < 0.0 ) def = 0.0
	782
	783	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	784
	785	ENDIF
	786
	787	ENDDO
	788	ENDDO
	789	ENDDO
	790
	791	IF ( prandtl_layer ) THEN
	792
	793	!
	794	!-- Position beneath wall
	795	!-- (2) - Will allways be executed.
	796	!-- 'bottom and wall: use u_0,v_0 and wall functions'
	797	!$acc loop
[1128]	798	DO i = i_left, i_right
	799	DO j = j_south, j_north
[1015]	800	!$acc loop vector( 32 )
	801	DO k = 1, nzt
	802
	803	IF ( ( wall_e_x(j,i) /= 0.0 ).OR.( wall_e_y(j,i) /= 0.0 ) ) &
	804	THEN
	805
	806	IF ( k == nzb_diff_s_inner(j,i) - 1 ) THEN
	807	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	808	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	809	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	810	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	811	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	812	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	813	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	814
	815	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	816	!
	817	!-- Inconsistency removed: as the thermal stratification is
	818	!-- not taken into account for the evaluation of the wall
	819	!-- fluxes at vertical walls, the eddy viscosity km must not
	820	!-- be used for the evaluation of the velocity gradients dudy
	821	!-- and dwdy
	822	!-- Note: The validity of the new method has not yet been
	823	!-- shown, as so far no suitable data for a validation
	824	!-- has been available
	825	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	826	! usvs, 1.0, 0.0, 0.0, 0.0 )
	827	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	828	! wsvs, 0.0, 0.0, 1.0, 0.0 )
	829	km_neutral = kappa * &
	830	( usvs(k,j,i)2 + wsvs(k,j,i)2 )*0.25 &
	831	0.5 * dy
	832	IF ( km_neutral > 0.0 ) THEN
	833	dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
	834	dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
	835	ELSE
	836	dudy = 0.0
	837	dwdy = 0.0
	838	ENDIF
	839	ELSE
	840	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	841	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	842	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	843	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	844	ENDIF
	845
	846	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	847	!
	848	!-- Inconsistency removed: as the thermal stratification is
	849	!-- not taken into account for the evaluation of the wall
	850	!-- fluxes at vertical walls, the eddy viscosity km must not
	851	!-- be used for the evaluation of the velocity gradients dvdx
	852	!-- and dwdx
	853	!-- Note: The validity of the new method has not yet been
	854	!-- shown, as so far no suitable data for a validation
	855	!-- has been available
	856	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	857	! vsus, 0.0, 1.0, 0.0, 0.0 )
	858	! CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	859	! wsus, 0.0, 0.0, 0.0, 1.0 )
	860	km_neutral = kappa * &
	861	( vsus(k,j,i)2 + wsus(k,j,i)2 )*0.25 &
	862	0.5 * dx
	863	IF ( km_neutral > 0.0 ) THEN
	864	dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
	865	dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
	866	ELSE
	867	dvdx = 0.0
	868	dwdx = 0.0
	869	ENDIF
	870	ELSE
	871	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	872	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	873	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	874	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	875	ENDIF
	876
	877	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	878	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	879	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	880
	881	IF ( def < 0.0 ) def = 0.0
	882
	883	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	884
	885	ENDIF
	886	!
	887	!-- (3) - will be executed only, if there is at least one level
	888	!-- between (2) and (4), i.e. the topography must have a
	889	!-- minimum height of 2 dz. Wall fluxes for this case have
	890	!-- already been calculated for (2).
	891	!-- 'wall only: use wall functions'
	892
	893	IF ( k >= nzb_diff_s_inner(j,i) .AND. &
	894	k <= nzb_diff_s_outer(j,i)-2 ) THEN
	895
	896	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	897	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	898	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	899	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	900	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	901	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	902	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	903
	904	IF ( wall_e_y(j,i) /= 0.0 ) THEN
	905	!
	906	!-- Inconsistency removed: as the thermal stratification
	907	!-- is not taken into account for the evaluation of the
	908	!-- wall fluxes at vertical walls, the eddy viscosity km
	909	!-- must not be used for the evaluation of the velocity
	910	!-- gradients dudy and dwdy
	911	!-- Note: The validity of the new method has not yet
	912	!-- been shown, as so far no suitable data for a
	913	!-- validation has been available
	914	km_neutral = kappa * ( usvs(k,j,i)**2 + &
	915	wsvs(k,j,i)2 )0.25 * 0.5 * dy
	916	IF ( km_neutral > 0.0 ) THEN
	917	dudy = - wall_e_y(j,i) * usvs(k,j,i) / km_neutral
	918	dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km_neutral
	919	ELSE
	920	dudy = 0.0
	921	dwdy = 0.0
	922	ENDIF
	923	ELSE
	924	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	925	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	926	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	927	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	928	ENDIF
	929
	930	IF ( wall_e_x(j,i) /= 0.0 ) THEN
	931	!
	932	!-- Inconsistency removed: as the thermal stratification
	933	!-- is not taken into account for the evaluation of the
	934	!-- wall fluxes at vertical walls, the eddy viscosity km
	935	!-- must not be used for the evaluation of the velocity
	936	!-- gradients dvdx and dwdx
	937	!-- Note: The validity of the new method has not yet
	938	!-- been shown, as so far no suitable data for a
	939	!-- validation has been available
	940	km_neutral = kappa * ( vsus(k,j,i)**2 + &
	941	wsus(k,j,i)2 )0.25 * 0.5 * dx
	942	IF ( km_neutral > 0.0 ) THEN
	943	dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km_neutral
	944	dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km_neutral
	945	ELSE
	946	dvdx = 0.0
	947	dwdx = 0.0
	948	ENDIF
	949	ELSE
	950	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	951	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	952	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	953	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	954	ENDIF
	955
	956	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	957	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	958	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	959
	960	IF ( def < 0.0 ) def = 0.0
	961
	962	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	963
	964	ENDIF
	965
	966	!
	967	!-- (4) - will allways be executed.
	968	!-- 'special case: free atmosphere' (as for case (0))
	969	IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN
	970
	971	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	972	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	973	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	974	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	975	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	976
	977	dvdx = 0.25 * ( 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	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	980	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	981	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	982
	983	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	984	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	985	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	986	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	987	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	988
	989	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	990	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	991	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	992
	993	IF ( def < 0.0 ) def = 0.0
	994
	995	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	996
	997	ENDIF
	998
	999	ENDIF
	1000
	1001	ENDDO
	1002	ENDDO
	1003	ENDDO
	1004
	1005	!
	1006	!-- Position without adjacent wall
	1007	!-- (1) - will allways be executed.
	1008	!-- 'bottom only: use u_0,v_0'
	1009	!$acc loop
[1128]	1010	DO i = i_left, i_right
	1011	DO j = j_south, j_north
[1015]	1012	!$acc loop vector( 32 )
	1013	DO k = 1, nzt
	1014
	1015	IF ( ( wall_e_x(j,i) == 0.0 ) .AND. ( wall_e_y(j,i) == 0.0 ) ) &
	1016	THEN
	1017
	1018	IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	1019
	1020	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1021	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1022	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1023	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1024	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	1025
	1026	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1027	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1028	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1029	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1030	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	1031
	1032	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1033	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1034	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1035	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1036	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1037
	1038	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	1039	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	1040	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1041
	1042	IF ( def < 0.0 ) def = 0.0
	1043
	1044	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1045
	1046	ENDIF
	1047
	1048	ENDIF
	1049
	1050	ENDDO
	1051	ENDDO
	1052	ENDDO
	1053
	1054	ELSEIF ( use_surface_fluxes ) THEN
	1055
	1056	!$acc loop
[1128]	1057	DO i = i_left, i_right
	1058	DO j = j_south, j_north
[1015]	1059	!$acc loop vector(32)
	1060	DO k = 1, nzt
	1061
	1062	IF ( k == nzb_diff_s_outer(j,i)-1 ) THEN
	1063
	1064	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1065	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1066	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1067	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1068	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	1069
	1070	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1071	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1072	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1073	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1074	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	1075
	1076	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1077	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1078	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1079	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1080	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1081
	1082	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	1083	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	1084	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1085
	1086	IF ( def < 0.0 ) def = 0.0
	1087
	1088	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1089
	1090	ENDIF
	1091
	1092	ENDDO
	1093	ENDDO
	1094	ENDDO
	1095
	1096	ENDIF
	1097
	1098	!
	1099	!-- If required, calculate TKE production by buoyancy
	1100	IF ( .NOT. neutral ) THEN
	1101
	1102	IF ( .NOT. humidity ) THEN
	1103
[1179]	1104	IF ( use_single_reference_value ) THEN
[1015]	1105
	1106	IF ( ocean ) THEN
	1107	!
	1108	!-- So far in the ocean no special treatment of density flux
	1109	!-- in the bottom and top surface layer
	1110	!$acc loop
[1128]	1111	DO i = i_left, i_right
	1112	DO j = j_south, j_north
[1015]	1113	!$acc loop vector( 32 )
	1114	DO k = 1, nzt
	1115	IF ( k > nzb_s_inner(j,i) ) THEN
	1116	tend(k,j,i) = tend(k,j,i) + &
	1117	kh(k,j,i) * g / rho_reference * &
	1118	( rho(k+1,j,i) - rho(k-1,j,i) ) * &
	1119	dd2zu(k)
	1120	ENDIF
	1121	ENDDO
	1122	ENDDO
	1123	ENDDO
	1124
	1125	ELSE
	1126
	1127	!$acc loop
[1128]	1128	DO i = i_left, i_right
	1129	DO j = j_south, j_north
[1015]	1130	!$acc loop vector( 32 )
	1131	DO k = 1, nzt_diff
	1132	IF ( k >= nzb_diff_s_inner(j,i) ) THEN
	1133	tend(k,j,i) = tend(k,j,i) - &
	1134	kh(k,j,i) * g / pt_reference * &
	1135	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
	1136	dd2zu(k)
	1137	ENDIF
	1138
	1139	IF ( k == nzb_diff_s_inner(j,i)-1 .AND. &
	1140	use_surface_fluxes ) THEN
	1141	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
	1142	shf(j,i)
	1143	ENDIF
	1144
	1145	IF ( k == nzt .AND. use_top_fluxes ) THEN
	1146	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
	1147	tswst(j,i)
	1148	ENDIF
	1149	ENDDO
	1150	ENDDO
	1151	ENDDO
	1152
	1153	ENDIF
	1154
	1155	ELSE
	1156
	1157	IF ( ocean ) THEN
	1158	!
	1159	!-- So far in the ocean no special treatment of density flux
	1160	!-- in the bottom and top surface layer
	1161	!$acc loop
[1128]	1162	DO i = i_left, i_right
	1163	DO j = j_south, j_north
[1015]	1164	!$acc loop vector( 32 )
	1165	DO k = 1, nzt
	1166	IF ( k > nzb_s_inner(j,i) ) THEN
	1167	tend(k,j,i) = tend(k,j,i) + &
	1168	kh(k,j,i) * g / rho(k,j,i) * &
	1169	( rho(k+1,j,i) - rho(k-1,j,i) ) * &
	1170	dd2zu(k)
	1171	ENDIF
	1172	ENDDO
	1173	ENDDO
	1174	ENDDO
	1175
	1176	ELSE
	1177
	1178	!$acc loop
[1128]	1179	DO i = i_left, i_right
	1180	DO j = j_south, j_north
[1015]	1181	!$acc loop vector( 32 )
	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(k,j,i) * &
	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(k,j,i) * &
	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(k,j,i) * &
	1198	tswst(j,i)
	1199	ENDIF
	1200	ENDDO
	1201	ENDDO
	1202	ENDDO
	1203
	1204	ENDIF
	1205
	1206	ENDIF
	1207
	1208	ELSE
	1209	!
	1210	!++ This part gives the PGI compiler problems in the previous loop
	1211	!++ even without any acc statements????
	1212	! STOP '+++ production_e problems with acc-directives'
	1213	! !acc loop
	1214	! DO i = nxl, nxr
	1215	! DO j = nys, nyn
	1216	! !acc loop vector( 32 )
	1217	! DO k = 1, nzt_diff
	1218	!
	1219	! IF ( k >= nzb_diff_s_inner(j,i) ) THEN
	1220	!
	1221	! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	1222	! k1 = 1.0 + 0.61 * q(k,j,i)
	1223	! k2 = 0.61 * pt(k,j,i)
	1224	! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	1225	! g / vpt(k,j,i) * &
	1226	! ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1227	! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1228	! ) * dd2zu(k)
	1229	! ELSE IF ( cloud_physics ) THEN
	1230	! IF ( ql(k,j,i) == 0.0 ) THEN
	1231	! k1 = 1.0 + 0.61 * q(k,j,i)
	1232	! k2 = 0.61 * pt(k,j,i)
	1233	! ELSE
	1234	! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1235	! temp = theta * t_d_pt(k)
	1236	! k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1237	! ( q(k,j,i) - ql(k,j,i) ) * &
	1238	! ( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1239	! ( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1240	! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1241	! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1242	! ENDIF
	1243	! tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	1244	! g / vpt(k,j,i) * &
	1245	! ( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1246	! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1247	! ) * dd2zu(k)
	1248	! ELSE IF ( cloud_droplets ) THEN
	1249	! k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	1250	! k2 = 0.61 * pt(k,j,i)
	1251	! tend(k,j,i) = tend(k,j,i) - &
	1252	! kh(k,j,i) * g / vpt(k,j,i) * &
	1253	! ( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &
	1254	! k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
	1255	! pt(k,j,i) * ( ql(k+1,j,i) - &
	1256	! ql(k-1,j,i) ) ) * dd2zu(k)
	1257	! ENDIF
	1258	!
	1259	! ENDIF
	1260	!
	1261	! ENDDO
	1262	! ENDDO
	1263	! ENDDO
	1264	!
	1265
	1266	!!++ Next two loops are probably very inefficiently parallellized
	1267	!!++ and will require better optimization
	1268	! IF ( use_surface_fluxes ) THEN
	1269	!
	1270	! !acc loop
	1271	! DO i = nxl, nxr
	1272	! DO j = nys, nyn
	1273	! !acc loop vector( 32 )
	1274	! DO k = 1, nzt_diff
	1275	!
	1276	! IF ( k == nzb_diff_s_inner(j,i)-1 ) THEN
	1277	!
	1278	! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	1279	! k1 = 1.0 + 0.61 * q(k,j,i)
	1280	! k2 = 0.61 * pt(k,j,i)
	1281	! ELSE IF ( cloud_physics ) THEN
	1282	! IF ( ql(k,j,i) == 0.0 ) THEN
	1283	! k1 = 1.0 + 0.61 * q(k,j,i)
	1284	! k2 = 0.61 * pt(k,j,i)
	1285	! ELSE
	1286	! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1287	! temp = theta * t_d_pt(k)
	1288	! k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1289	! ( q(k,j,i) - ql(k,j,i) ) * &
	1290	! ( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1291	! ( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1292	! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1293	! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1294	! ENDIF
	1295	! ELSE IF ( cloud_droplets ) THEN
	1296	! k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	1297	! k2 = 0.61 * pt(k,j,i)
	1298	! ENDIF
	1299	!
	1300	! tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1301	! ( k1* shf(j,i) + k2 * qsws(j,i) )
	1302	! ENDIF
	1303	!
	1304	! ENDDO
	1305	! ENDDO
	1306	! ENDDO
	1307	!
	1308	! ENDIF
	1309	!
	1310	! IF ( use_top_fluxes ) THEN
	1311	!
	1312	! !acc loop
	1313	! DO i = nxl, nxr
	1314	! DO j = nys, nyn
	1315	! !acc loop vector( 32 )
	1316	! DO k = 1, nzt
	1317	! IF ( k == nzt ) THEN
	1318	!
	1319	! IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
	1320	! k1 = 1.0 + 0.61 * q(k,j,i)
	1321	! k2 = 0.61 * pt(k,j,i)
	1322	! ELSE IF ( cloud_physics ) THEN
	1323	! IF ( ql(k,j,i) == 0.0 ) THEN
	1324	! k1 = 1.0 + 0.61 * q(k,j,i)
	1325	! k2 = 0.61 * pt(k,j,i)
	1326	! ELSE
	1327	! theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1328	! temp = theta * t_d_pt(k)
	1329	! k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1330	! ( q(k,j,i) - ql(k,j,i) ) * &
	1331	! ( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1332	! ( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1333	! ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1334	! k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1335	! ENDIF
	1336	! ELSE IF ( cloud_droplets ) THEN
	1337	! k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	1338	! k2 = 0.61 * pt(k,j,i)
	1339	! ENDIF
	1340	!
	1341	! tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1342	! ( k1* tswst(j,i) + k2 * qswst(j,i) )
	1343	!
	1344	! ENDIF
	1345	!
	1346	! ENDDO
	1347	! ENDDO
	1348	! ENDDO
	1349	!
	1350	! ENDIF
	1351
	1352	ENDIF
	1353
	1354	ENDIF
	1355	!$acc end kernels
	1356
	1357	END SUBROUTINE production_e_acc
	1358
	1359
	1360	!------------------------------------------------------------------------------!
[1]	1361	! Call for grid point i,j
	1362	!------------------------------------------------------------------------------!
	1363	SUBROUTINE production_e_ij( i, j )
	1364
	1365	USE arrays_3d
	1366	USE cloud_parameters
	1367	USE control_parameters
	1368	USE grid_variables
	1369	USE indices
	1370	USE statistics
[449]	1371
[1]	1372	IMPLICIT NONE
	1373
	1374	INTEGER :: i, j, k
	1375
	1376	REAL :: def, dudx, dudy, dudz, dvdx, dvdy, dvdz, dwdx, dwdy, dwdz, &
[208]	1377	k1, k2, km_neutral, theta, temp
[1]	1378
[56]	1379	REAL, DIMENSION(nzb:nzt+1) :: usvs, vsus, wsus, wsvs
[53]	1380
[1]	1381	!
	1382	!-- Calculate TKE production by shear
[19]	1383	DO k = nzb_diff_s_outer(j,i), nzt
[1]	1384
	1385	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1386	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1387	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1388	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1389	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	1390
	1391	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1392	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1393	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1394	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1395	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	1396
	1397	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1398	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1399	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1400	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1401	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1402
	1403	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	1404	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	1405	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1406
	1407	IF ( def < 0.0 ) def = 0.0
	1408
	1409	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
[1007]	1410
[1]	1411	ENDDO
	1412
[37]	1413	IF ( prandtl_layer ) THEN
[1]	1414
	1415	IF ( ( wall_e_x(j,i) /= 0.0 ) .OR. ( wall_e_y(j,i) /= 0.0 ) ) THEN
[55]	1416
[1]	1417	!
[55]	1418	!-- Position beneath wall
	1419	!-- (2) - Will allways be executed.
	1420	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	1421	k = nzb_diff_s_inner(j,i)-1
	1422
	1423	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	1424	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1425	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	1426	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1427	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1428	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	1429	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1430
[1]	1431	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[1007]	1432	!
[208]	1433	!-- Inconsistency removed: as the thermal stratification
	1434	!-- is not taken into account for the evaluation of the
	1435	!-- wall fluxes at vertical walls, the eddy viscosity km
	1436	!-- must not be used for the evaluation of the velocity
	1437	!-- gradients dudy and dwdy
	1438	!-- Note: The validity of the new method has not yet
	1439	!-- been shown, as so far no suitable data for a
	1440	!-- validation has been available
[53]	1441	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	1442	usvs, 1.0, 0.0, 0.0, 0.0 )
	1443	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	1444	wsvs, 0.0, 0.0, 1.0, 0.0 )
[208]	1445	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25 &
	1446	0.5 * dy
[364]	1447	IF ( km_neutral > 0.0 ) THEN
	1448	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	1449	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	1450	ELSE
	1451	dudy = 0.0
	1452	dwdy = 0.0
	1453	ENDIF
[1]	1454	ELSE
	1455	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1456	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	1457	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1458	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	1459	ENDIF
	1460
	1461	IF ( wall_e_x(j,i) /= 0.0 ) THEN
[1007]	1462	!
[208]	1463	!-- Inconsistency removed: as the thermal stratification
	1464	!-- is not taken into account for the evaluation of the
	1465	!-- wall fluxes at vertical walls, the eddy viscosity km
	1466	!-- must not be used for the evaluation of the velocity
	1467	!-- gradients dvdx and dwdx
	1468	!-- Note: The validity of the new method has not yet
	1469	!-- been shown, as so far no suitable data for a
	1470	!-- validation has been available
[53]	1471	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	1472	vsus, 0.0, 1.0, 0.0, 0.0 )
	1473	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
	1474	wsus, 0.0, 0.0, 0.0, 1.0 )
[208]	1475	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25 &
	1476	0.5 * dx
[364]	1477	IF ( km_neutral > 0.0 ) THEN
	1478	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	1479	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	1480	ELSE
	1481	dvdx = 0.0
	1482	dwdx = 0.0
	1483	ENDIF
[1]	1484	ELSE
	1485	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1486	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1487	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1488	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1489	ENDIF
	1490
	1491	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	1492	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	1493	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1494
	1495	IF ( def < 0.0 ) def = 0.0
	1496
	1497	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1498
	1499	!
[55]	1500	!-- (3) - will be executed only, if there is at least one level
	1501	!-- between (2) and (4), i.e. the topography must have a
	1502	!-- minimum height of 2 dz. Wall fluxes for this case have
	1503	!-- already been calculated for (2).
	1504	!-- 'wall only: use wall functions'
[1]	1505	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	1506
	1507	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[53]	1508	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1509	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	1510	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1511	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1512	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	1513	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1514
[1]	1515	IF ( wall_e_y(j,i) /= 0.0 ) THEN
[1007]	1516	!
[208]	1517	!-- Inconsistency removed: as the thermal stratification
	1518	!-- is not taken into account for the evaluation of the
	1519	!-- wall fluxes at vertical walls, the eddy viscosity km
	1520	!-- must not be used for the evaluation of the velocity
	1521	!-- gradients dudy and dwdy
	1522	!-- Note: The validity of the new method has not yet
	1523	!-- been shown, as so far no suitable data for a
	1524	!-- validation has been available
	1525	km_neutral = kappa * ( usvs(k)**2 + &
	1526	wsvs(k)2 )0.25 * 0.5 * dy
[364]	1527	IF ( km_neutral > 0.0 ) THEN
	1528	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	1529	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	1530	ELSE
	1531	dudy = 0.0
	1532	dwdy = 0.0
	1533	ENDIF
[1]	1534	ELSE
	1535	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1536	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
[53]	1537	dwdy = 0.25 * ( 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
	1541	IF ( wall_e_x(j,i) /= 0.0 ) 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
	1551	km_neutral = kappa * ( vsus(k)**2 + &
	1552	wsus(k)2 )0.25 * 0.5 * dx
[364]	1553	IF ( km_neutral > 0.0 ) THEN
	1554	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	1555	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	1556	ELSE
	1557	dvdx = 0.0
	1558	dwdx = 0.0
	1559	ENDIF
[1]	1560	ELSE
	1561	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1562	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1563	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1564	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1565	ENDIF
	1566
	1567	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	1568	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	1569	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1570
	1571	IF ( def < 0.0 ) def = 0.0
	1572
	1573	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1574
	1575	ENDDO
	1576
	1577	!
[55]	1578	!-- (4) - will allways be executed.
	1579	!-- 'special case: free atmosphere' (as for case (0))
[1]	1580	k = nzb_diff_s_outer(j,i)-1
	1581
	1582	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1583	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1584	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1585	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1586	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	1587
	1588	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1589	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1590	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1591	dvdz = 0.5 * ( 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)
	1593
	1594	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1595	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1596	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1597	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1598	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1599
	1600	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	1601	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	1602	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1603
	1604	IF ( def < 0.0 ) def = 0.0
	1605
	1606	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1607
	1608	ELSE
	1609
	1610	!
[55]	1611	!-- Position without adjacent wall
	1612	!-- (1) - will allways be executed.
	1613	!-- 'bottom only: use u_0,v_0'
[1]	1614	k = nzb_diff_s_inner(j,i)-1
	1615
	1616	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1617	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1618	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1619	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1620	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	1621
	1622	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1623	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1624	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1625	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1626	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
	1627
	1628	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1629	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1630	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1631	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1632	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1633
	1634	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) &
	1635	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	1636	+ 2.0 * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1637
	1638	IF ( def < 0.0 ) def = 0.0
	1639
	1640	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1641
	1642	ENDIF
	1643
[37]	1644	ELSEIF ( use_surface_fluxes ) THEN
	1645
	1646	k = nzb_diff_s_outer(j,i)-1
	1647
	1648	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1649	dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1650	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1651	dudz = 0.5 * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1652	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	1653
	1654	dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1655	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1656	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1657	dvdz = 0.5 * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1658	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
	1659
	1660	dwdx = 0.25 * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1661	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1662	dwdy = 0.25 * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1663	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1664	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	1665
	1666	def = 2.0 * ( dudx2 + dvdy2 + dwdz**2 ) + &
	1667	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
	1668	dvdz*2 + 2.0 ( dvdxdudy + dwdxdudz + dwdy*dvdz )
	1669
	1670	IF ( def < 0.0 ) def = 0.0
	1671
	1672	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1673
[1]	1674	ENDIF
	1675
	1676	!
[940]	1677	!-- If required, calculate TKE production by buoyancy
	1678	IF ( .NOT. neutral ) THEN
[1]	1679
[940]	1680	IF ( .NOT. humidity ) THEN
[19]	1681
[1179]	1682	IF ( use_single_reference_value ) THEN
[940]	1683
	1684	IF ( ocean ) THEN
[97]	1685	!
[940]	1686	!-- So far in the ocean no special treatment of density flux in
	1687	!-- the bottom and top surface layer
	1688	DO k = nzb_s_inner(j,i)+1, nzt
	1689	tend(k,j,i) = tend(k,j,i) + &
	1690	kh(k,j,i) * g / rho_reference * &
	1691	( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
	1692	ENDDO
[97]	1693
[940]	1694	ELSE
[97]	1695
[940]	1696	DO k = nzb_diff_s_inner(j,i), nzt_diff
	1697	tend(k,j,i) = tend(k,j,i) - &
	1698	kh(k,j,i) * g / pt_reference * &
	1699	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	1700	ENDDO
[1]	1701
[940]	1702	IF ( use_surface_fluxes ) THEN
	1703	k = nzb_diff_s_inner(j,i)-1
	1704	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
	1705	ENDIF
[19]	1706
[940]	1707	IF ( use_top_fluxes ) THEN
	1708	k = nzt
	1709	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
	1710	ENDIF
	1711
[97]	1712	ENDIF
	1713
[940]	1714	ELSE
[57]	1715
[940]	1716	IF ( ocean ) THEN
[97]	1717	!
[940]	1718	!-- So far in the ocean no special treatment of density flux in
	1719	!-- the bottom and top surface layer
	1720	DO k = nzb_s_inner(j,i)+1, nzt
	1721	tend(k,j,i) = tend(k,j,i) + &
	1722	kh(k,j,i) * g / rho(k,j,i) * &
	1723	( rho(k+1,j,i) - rho(k-1,j,i) ) * dd2zu(k)
	1724	ENDDO
[97]	1725
[940]	1726	ELSE
[97]	1727
[940]	1728	DO k = nzb_diff_s_inner(j,i), nzt_diff
	1729	tend(k,j,i) = tend(k,j,i) - &
	1730	kh(k,j,i) * g / pt(k,j,i) * &
	1731	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	1732	ENDDO
[57]	1733
[940]	1734	IF ( use_surface_fluxes ) THEN
	1735	k = nzb_diff_s_inner(j,i)-1
	1736	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	1737	ENDIF
[57]	1738
[940]	1739	IF ( use_top_fluxes ) THEN
	1740	k = nzt
	1741	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	1742	ENDIF
	1743
[97]	1744	ENDIF
	1745
[57]	1746	ENDIF
	1747
[940]	1748	ELSE
[57]	1749
[940]	1750	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	1751
[1007]	1752	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1]	1753	k1 = 1.0 + 0.61 * q(k,j,i)
	1754	k2 = 0.61 * pt(k,j,i)
[1007]	1755	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1756	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1757	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1758	) * dd2zu(k)
	1759	ELSE IF ( cloud_physics ) THEN
[940]	1760	IF ( ql(k,j,i) == 0.0 ) THEN
	1761	k1 = 1.0 + 0.61 * q(k,j,i)
	1762	k2 = 0.61 * pt(k,j,i)
	1763	ELSE
	1764	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1765	temp = theta * t_d_pt(k)
	1766	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1767	( q(k,j,i) - ql(k,j,i) ) * &
	1768	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1769	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1770	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1771	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1772	ENDIF
[1007]	1773	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
[940]	1774	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1775	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1776	) * dd2zu(k)
[1007]	1777	ELSE IF ( cloud_droplets ) THEN
	1778	k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	1779	k2 = 0.61 * pt(k,j,i)
	1780	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1781	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1782	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
	1783	pt(k,j,i) * ( ql(k+1,j,i) - &
	1784	ql(k-1,j,i) ) ) * dd2zu(k)
	1785	ENDIF
[940]	1786	ENDDO
[19]	1787
[940]	1788	IF ( use_surface_fluxes ) THEN
	1789	k = nzb_diff_s_inner(j,i)-1
[1]	1790
[1007]	1791	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1]	1792	k1 = 1.0 + 0.61 * q(k,j,i)
	1793	k2 = 0.61 * pt(k,j,i)
[1007]	1794	ELSE IF ( cloud_physics ) THEN
[940]	1795	IF ( ql(k,j,i) == 0.0 ) THEN
	1796	k1 = 1.0 + 0.61 * q(k,j,i)
	1797	k2 = 0.61 * pt(k,j,i)
	1798	ELSE
	1799	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1800	temp = theta * t_d_pt(k)
	1801	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1802	( q(k,j,i) - ql(k,j,i) ) * &
	1803	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1804	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1805	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1806	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1807	ENDIF
[1007]	1808	ELSE IF ( cloud_droplets ) THEN
	1809	k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	1810	k2 = 0.61 * pt(k,j,i)
[1]	1811	ENDIF
[940]	1812
	1813	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1814	( k1* shf(j,i) + k2 * qsws(j,i) )
[1]	1815	ENDIF
	1816
[940]	1817	IF ( use_top_fluxes ) THEN
	1818	k = nzt
[1]	1819
[1007]	1820	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[19]	1821	k1 = 1.0 + 0.61 * q(k,j,i)
	1822	k2 = 0.61 * pt(k,j,i)
[1007]	1823	ELSE IF ( cloud_physics ) THEN
[940]	1824	IF ( ql(k,j,i) == 0.0 ) THEN
	1825	k1 = 1.0 + 0.61 * q(k,j,i)
	1826	k2 = 0.61 * pt(k,j,i)
	1827	ELSE
	1828	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1829	temp = theta * t_d_pt(k)
	1830	k1 = ( 1.0 - q(k,j,i) + 1.61 * &
	1831	( q(k,j,i) - ql(k,j,i) ) * &
	1832	( 1.0 + 0.622 * l_d_r / temp ) ) / &
	1833	( 1.0 + 0.622 * l_d_r * l_d_cp * &
	1834	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
	1835	k2 = theta * ( l_d_cp / temp * k1 - 1.0 )
	1836	ENDIF
[1007]	1837	ELSE IF ( cloud_droplets ) THEN
	1838	k1 = 1.0 + 0.61 * q(k,j,i) - ql(k,j,i)
	1839	k2 = 0.61 * pt(k,j,i)
[19]	1840	ENDIF
[940]	1841
	1842	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1843	( k1* tswst(j,i) + k2 * qswst(j,i) )
[19]	1844	ENDIF
	1845
	1846	ENDIF
	1847
[1]	1848	ENDIF
	1849
	1850	END SUBROUTINE production_e_ij
	1851
	1852
	1853	SUBROUTINE production_e_init
	1854
	1855	USE arrays_3d
	1856	USE control_parameters
	1857	USE grid_variables
	1858	USE indices
	1859
	1860	IMPLICIT NONE
	1861
	1862	INTEGER :: i, j, ku, kv
	1863
[37]	1864	IF ( prandtl_layer ) THEN
[1]	1865
	1866	IF ( first_call ) THEN
[759]	1867	ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
	1868	u_0 = 0.0 ! just to avoid access of uninitialized memory
	1869	v_0 = 0.0 ! within exchange_horiz_2d
[1]	1870	first_call = .FALSE.
	1871	ENDIF
	1872
	1873	!
	1874	!-- Calculate a virtual velocity at the surface in a way that the
	1875	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
	1876	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
	1877	!-- production term at k=1 (see production_e_ij).
	1878	!-- The velocity gradient has to be limited in case of too small km
	1879	!-- (otherwise the timestep may be significantly reduced by large
	1880	!-- surface winds).
[106]	1881	!-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are
	1882	!-- not available in case of non-cyclic boundary conditions.
[1]	1883	!-- WARNING: the exact analytical solution would require the determination
	1884	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
	1885	!$OMP PARALLEL DO PRIVATE( ku, kv )
[106]	1886	DO i = nxl, nxr+1
	1887	DO j = nys, nyn+1
[1]	1888
	1889	ku = nzb_u_inner(j,i)+1
	1890	kv = nzb_v_inner(j,i)+1
	1891
	1892	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
	1893	( 0.5 * ( km(ku,j,i) + km(ku,j,i-1) ) + &
	1894	1.0E-20 )
	1895	! ( us(j,i) * kappa * zu(1) )
	1896	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
	1897	( 0.5 * ( km(kv,j,i) + km(kv,j-1,i) ) + &
	1898	1.0E-20 )
	1899	! ( us(j,i) * kappa * zu(1) )
	1900
	1901	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
	1902	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
	1903	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
	1904	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
	1905
	1906	ENDDO
	1907	ENDDO
	1908
	1909	CALL exchange_horiz_2d( u_0 )
	1910	CALL exchange_horiz_2d( v_0 )
	1911
	1912	ENDIF
	1913
	1914	END SUBROUTINE production_e_init
	1915
	1916	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |