Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |