Home

Context Navigation

production_e.f90 @ 1180

Last change on this file since 1180 was 1179, checked in by raasch, 12 years ago

New:
---
Initial profiles can be used as reference state in the buoyancy term. New parameter
reference_state introduced. Calculation and handling of reference state in buoyancy term revised.
binary version for restart files changed from 3.9 to 3.9a (no downward compatibility!),
initial profile for rho added to hom (id=77)

Errors:

small bugfix for background communication (time_integration)

Property svn:keywords set to Id

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

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

Context Navigation

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

Download in other formats: