Home

Context Navigation

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

Last change on this file since 1267 was 1258, checked in by raasch, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 82.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |