Home

Context Navigation

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

Last change on this file since 1684 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 89.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |