Home

Context Navigation

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

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

Rev	Line
[1873]	1	!> @file production_e.f90
[2000]	2	!------------------------------------------------------------------------------!
[1036]	3	! This file is part of PALM.
	4	!
[2000]	5	! PALM is free software: you can redistribute it and/or modify it under the
	6	! terms of the GNU General Public License as published by the Free Software
	7	! Foundation, either version 3 of the License, or (at your option) any later
	8	! version.
[1036]	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	!
[2101]	17	! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]	18	!------------------------------------------------------------------------------!
[1036]	19	!
[484]	20	! Current revisions:
[1]	21	! -----------------
[1343]	22	!
[2119]	23	!
[1321]	24	! Former revisions:
	25	! -----------------
	26	! $Id: production_e.f90 2119 2017-01-17 16:51:50Z raasch $
	27	!
[2119]	28	! 2118 2017-01-17 16:38:49Z raasch
	29	! OpenACC version of subroutine removed
	30	!
[2032]	31	! 2031 2016-10-21 15:11:58Z knoop
	32	! renamed variable rho to rho_ocean
	33	!
[2001]	34	! 2000 2016-08-20 18:09:15Z knoop
	35	! Forced header and separation lines into 80 columns
	36	!
[1874]	37	! 1873 2016-04-18 14:50:06Z maronga
	38	! Module renamed (removed _mod)
	39	!
	40	!
[1851]	41	! 1850 2016-04-08 13:29:27Z maronga
	42	! Module renamed
	43	!
	44	!
[1692]	45	! 1691 2015-10-26 16:17:44Z maronga
	46	! Renamed prandtl_layer to constant_flux_layer.
	47	!
[1683]	48	! 1682 2015-10-07 23:56:08Z knoop
	49	! Code annotations made doxygen readable
	50	!
[1375]	51	! 1374 2014-04-25 12:55:07Z raasch
	52	! nzb_s_outer removed from acc-present-list
	53	!
[1354]	54	! 1353 2014-04-08 15:21:23Z heinze
	55	! REAL constants provided with KIND-attribute
	56	!
[1343]	57	! 1342 2014-03-26 17:04:47Z kanani
	58	! REAL constants defined as wp-kind
	59	!
[1321]	60	! 1320 2014-03-20 08:40:49Z raasch
[1320]	61	! ONLY-attribute added to USE-statements,
	62	! kind-parameters added to all INTEGER and REAL declaration statements,
	63	! kinds are defined in new module kinds,
	64	! old module precision_kind is removed,
	65	! revision history before 2012 removed,
	66	! comment fields (!:) to be used for variable explanations added to
	67	! all variable declaration statements
[110]	68	!
[1258]	69	! 1257 2013-11-08 15:18:40Z raasch
	70	! openacc loop and loop vector clauses removed, declare create moved after
	71	! the FORTRAN declaration statement
	72	!
[1182]	73	! 1179 2013-06-14 05:57:58Z raasch
	74	! use_reference renamed use_single_reference_value
	75	!
[1132]	76	! 1128 2013-04-12 06:19:32Z raasch
	77	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	78	! j_north
	79	!
[1037]	80	! 1036 2012-10-22 13:43:42Z raasch
	81	! code put under GPL (PALM 3.9)
	82	!
[1017]	83	! 1015 2012-09-27 09:23:24Z raasch
	84	! accelerator version (*_acc) added
	85	!
[1008]	86	! 1007 2012-09-19 14:30:36Z franke
	87	! Bugfix: calculation of buoyancy production has to consider the liquid water
	88	! mixing ratio in case of cloud droplets
	89	!
[941]	90	! 940 2012-07-09 14:31:00Z raasch
	91	! TKE production by buoyancy can be switched off in case of runs with pure
	92	! neutral stratification
	93	!
[1]	94	! Revision 1.1 1997/09/19 07:45:35 raasch
	95	! Initial revision
	96	!
	97	!
	98	! Description:
	99	! ------------
[1682]	100	!> Production terms (shear + buoyancy) of the TKE.
[1691]	101	!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
[1682]	102	!> not considered well!
[1]	103	!------------------------------------------------------------------------------!
[1682]	104	MODULE production_e_mod
	105
[1]	106
[1320]	107	USE wall_fluxes_mod, &
[2118]	108	ONLY: wall_fluxes_e
[56]	109
[1320]	110	USE kinds
	111
[1]	112	PRIVATE
[2118]	113	PUBLIC production_e, production_e_init
[56]	114
[1682]	115	LOGICAL, SAVE :: first_call = .TRUE. !<
[1]	116
[1682]	117	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: u_0 !<
	118	REAL(wp), DIMENSION(:,:), ALLOCATABLE, SAVE :: v_0 !<
[1]	119
	120	INTERFACE production_e
	121	MODULE PROCEDURE production_e
	122	MODULE PROCEDURE production_e_ij
	123	END INTERFACE production_e
	124
	125	INTERFACE production_e_init
	126	MODULE PROCEDURE production_e_init
	127	END INTERFACE production_e_init
	128
	129	CONTAINS
	130
	131
	132	!------------------------------------------------------------------------------!
[1682]	133	! Description:
	134	! ------------
	135	!> Call for all grid points
[1]	136	!------------------------------------------------------------------------------!
	137	SUBROUTINE production_e
	138
[1320]	139	USE arrays_3d, &
[2031]	140	ONLY: ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho_ocean, shf, &
[1320]	141	tend, tswst, u, v, vpt, w
[1]	142
[1320]	143	USE cloud_parameters, &
	144	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
	145
	146	USE control_parameters, &
[1691]	147	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
	148	humidity, kappa, neutral, ocean, pt_reference, &
	149	rho_reference, use_single_reference_value, &
	150	use_surface_fluxes, use_top_fluxes
[1320]	151
	152	USE grid_variables, &
	153	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
	154
	155	USE indices, &
	156	ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, &
	157	nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
	158
[1]	159	IMPLICIT NONE
	160
[1682]	161	INTEGER(iwp) :: i !<
	162	INTEGER(iwp) :: j !<
	163	INTEGER(iwp) :: k !<
[1]	164
[1682]	165	REAL(wp) :: def !<
	166	REAL(wp) :: dudx !<
	167	REAL(wp) :: dudy !<
	168	REAL(wp) :: dudz !<
	169	REAL(wp) :: dvdx !<
	170	REAL(wp) :: dvdy !<
	171	REAL(wp) :: dvdz !<
	172	REAL(wp) :: dwdx !<
	173	REAL(wp) :: dwdy !<
	174	REAL(wp) :: dwdz !<
	175	REAL(wp) :: k1 !<
	176	REAL(wp) :: k2 !<
	177	REAL(wp) :: km_neutral !<
	178	REAL(wp) :: theta !<
	179	REAL(wp) :: temp !<
[1]	180
[1320]	181	! REAL(wp), DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) :: usvs, vsus, wsus, wsvs
[1682]	182	REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !<
	183	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
	184	REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !<
	185	REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !<
[1]	186
[56]	187	!
	188	!-- First calculate horizontal momentum flux u'v', w'v', v'u', w'u' at
	189	!-- vertical walls, if neccessary
	190	!-- So far, results are slightly different from the ij-Version.
	191	!-- Therefore, ij-Version is called further below within the ij-loops.
	192	! IF ( topography /= 'flat' ) THEN
[1320]	193	! CALL wall_fluxes_e( usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, wall_e_y )
	194	! CALL wall_fluxes_e( wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp, wall_e_y )
	195	! CALL wall_fluxes_e( vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, wall_e_x )
	196	! CALL wall_fluxes_e( wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp, wall_e_x )
[56]	197	! ENDIF
[53]	198
[940]	199
[1]	200	DO i = nxl, nxr
	201
[940]	202	!
	203	!-- Calculate TKE production by shear
[1]	204	DO j = nys, nyn
[19]	205	DO k = nzb_diff_s_outer(j,i), nzt
[1]	206
[1342]	207	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	208	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	209	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	210	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	211	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[1]	212
[1342]	213	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	214	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	215	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	216	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	217	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[1]	218
[1342]	219	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	220	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	221	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	222	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	223	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	224
[1342]	225	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	226	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	227	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	228
[1342]	229	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	230
	231	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
[1007]	232
[1]	233	ENDDO
	234	ENDDO
	235
[1691]	236	IF ( constant_flux_layer ) THEN
[1]	237
	238	!
[55]	239	!-- Position beneath wall
	240	!-- (2) - Will allways be executed.
	241	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	242	DO j = nys, nyn
	243
[1342]	244	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
[1]	245	THEN
	246
	247	k = nzb_diff_s_inner(j,i) - 1
	248	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[1342]	249	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	250	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
[53]	251	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
[1342]	252	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	253	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
[53]	254	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	255
[1342]	256	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
[1007]	257	!
[208]	258	!-- Inconsistency removed: as the thermal stratification is
	259	!-- not taken into account for the evaluation of the wall
	260	!-- fluxes at vertical walls, the eddy viscosity km must not
	261	!-- be used for the evaluation of the velocity gradients dudy
	262	!-- and dwdy
	263	!-- Note: The validity of the new method has not yet been
	264	!-- shown, as so far no suitable data for a validation
	265	!-- has been available
[53]	266	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	267	usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
[53]	268	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	269	wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
[1342]	270	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25_wp &
	271	0.5_wp * dy
	272	IF ( km_neutral > 0.0_wp ) THEN
[364]	273	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	274	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	275	ELSE
[1342]	276	dudy = 0.0_wp
	277	dwdy = 0.0_wp
[364]	278	ENDIF
[1]	279	ELSE
[1342]	280	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	281	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	282	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	283	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	284	ENDIF
	285
[1342]	286	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
[1007]	287	!
[208]	288	!-- Inconsistency removed: as the thermal stratification is
	289	!-- not taken into account for the evaluation of the wall
	290	!-- fluxes at vertical walls, the eddy viscosity km must not
	291	!-- be used for the evaluation of the velocity gradients dvdx
	292	!-- and dwdx
	293	!-- Note: The validity of the new method has not yet been
	294	!-- shown, as so far no suitable data for a validation
	295	!-- has been available
[53]	296	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	297	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
[53]	298	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	299	wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
[1342]	300	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25_wp &
	301	0.5_wp * dx
	302	IF ( km_neutral > 0.0_wp ) THEN
[364]	303	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	304	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	305	ELSE
[1342]	306	dvdx = 0.0_wp
	307	dwdx = 0.0_wp
[364]	308	ENDIF
[1]	309	ELSE
[1342]	310	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	311	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	312	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	313	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
[1]	314	ENDIF
	315
[1342]	316	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	317	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	318	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	319
[1342]	320	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	321
	322	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	323
	324
	325	!
[55]	326	!-- (3) - will be executed only, if there is at least one level
	327	!-- between (2) and (4), i.e. the topography must have a
	328	!-- minimum height of 2 dz. Wall fluxes for this case have
	329	!-- already been calculated for (2).
	330	!-- 'wall only: use wall functions'
[1]	331
	332	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	333
	334	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[1342]	335	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	336	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	337	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	338	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	339	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[53]	340	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	341
[1342]	342	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
[1007]	343	!
[208]	344	!-- Inconsistency removed: as the thermal stratification
	345	!-- is not taken into account for the evaluation of the
	346	!-- wall fluxes at vertical walls, the eddy viscosity km
	347	!-- must not be used for the evaluation of the velocity
	348	!-- gradients dudy and dwdy
	349	!-- Note: The validity of the new method has not yet
	350	!-- been shown, as so far no suitable data for a
	351	!-- validation has been available
	352	km_neutral = kappa * ( usvs(k)**2 + &
[1342]	353	wsvs(k)2 )0.25_wp * 0.5_wp * dy
	354	IF ( km_neutral > 0.0_wp ) THEN
[364]	355	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	356	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	357	ELSE
[1342]	358	dudy = 0.0_wp
	359	dwdy = 0.0_wp
[364]	360	ENDIF
[1]	361	ELSE
[1342]	362	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	363	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	364	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	365	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	366	ENDIF
	367
[1342]	368	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
[1007]	369	!
[208]	370	!-- Inconsistency removed: as the thermal stratification
	371	!-- is not taken into account for the evaluation of the
	372	!-- wall fluxes at vertical walls, the eddy viscosity km
	373	!-- must not be used for the evaluation of the velocity
	374	!-- gradients dvdx and dwdx
	375	!-- Note: The validity of the new method has not yet
	376	!-- been shown, as so far no suitable data for a
	377	!-- validation has been available
	378	km_neutral = kappa * ( vsus(k)**2 + &
[1342]	379	wsus(k)2 )0.25_wp * 0.5_wp * dx
	380	IF ( km_neutral > 0.0_wp ) THEN
[364]	381	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	382	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	383	ELSE
[1342]	384	dvdx = 0.0_wp
	385	dwdx = 0.0_wp
[364]	386	ENDIF
[1]	387	ELSE
[1342]	388	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	389	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	390	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	391	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
[1]	392	ENDIF
	393
[1342]	394	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	395	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	396	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	397
[1342]	398	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	399
	400	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	401
	402	ENDDO
	403
	404	ENDIF
	405
	406	ENDDO
	407
	408	!
[55]	409	!-- (4) - will allways be executed.
	410	!-- 'special case: free atmosphere' (as for case (0))
[1]	411	DO j = nys, nyn
	412
[1342]	413	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) &
[1]	414	THEN
	415
	416	k = nzb_diff_s_outer(j,i)-1
	417
[1342]	418	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	419	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	420	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	421	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	422	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[1]	423
[1342]	424	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	425	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	426	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	427	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	428	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[1]	429
[1342]	430	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	431	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	432	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	433	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	434	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	435
[1342]	436	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	437	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	438	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	439
[1342]	440	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	441
	442	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	443
	444	ENDIF
	445
	446	ENDDO
	447
	448	!
[55]	449	!-- Position without adjacent wall
	450	!-- (1) - will allways be executed.
	451	!-- 'bottom only: use u_0,v_0'
[1]	452	DO j = nys, nyn
	453
[1342]	454	IF ( ( wall_e_x(j,i) == 0.0_wp ) .AND. ( wall_e_y(j,i) == 0.0_wp ) ) &
[1]	455	THEN
	456
	457	k = nzb_diff_s_inner(j,i)-1
	458
[1342]	459	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	460	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	461	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	462	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	463	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
[1]	464
[1342]	465	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	466	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	467	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	468	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	469	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
[1]	470
[1342]	471	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	472	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	473	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	474	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	475	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	476
[1342]	477	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	478	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	479	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	480
[1342]	481	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	482
	483	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
[1007]	484
[1]	485	ENDIF
	486
	487	ENDDO
	488
[37]	489	ELSEIF ( use_surface_fluxes ) THEN
	490
	491	DO j = nys, nyn
	492
	493	k = nzb_diff_s_outer(j,i)-1
	494
[1342]	495	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	496	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	497	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	498	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	499	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[37]	500
[1342]	501	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	502	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	503	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	504	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	505	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[37]	506
[1342]	507	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	508	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	509	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	510	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	511	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[37]	512
[1342]	513	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[37]	514	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	515	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[37]	516
[1342]	517	IF ( def < 0.0_wp ) def = 0.0_wp
[37]	518
	519	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	520
	521	ENDDO
	522
[1]	523	ENDIF
	524
	525	!
[940]	526	!-- If required, calculate TKE production by buoyancy
	527	IF ( .NOT. neutral ) THEN
[1]	528
[940]	529	IF ( .NOT. humidity ) THEN
[1]	530
[1179]	531	IF ( use_single_reference_value ) THEN
[940]	532
	533	IF ( ocean ) THEN
[97]	534	!
[940]	535	!-- So far in the ocean no special treatment of density flux
	536	!-- in the bottom and top surface layer
	537	DO j = nys, nyn
	538	DO k = nzb_s_inner(j,i)+1, nzt
	539	tend(k,j,i) = tend(k,j,i) + &
	540	kh(k,j,i) * g / rho_reference * &
[2031]	541	( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &
[940]	542	dd2zu(k)
	543	ENDDO
[97]	544	ENDDO
	545
[940]	546	ELSE
[97]	547
[940]	548	DO j = nys, nyn
	549	DO k = nzb_diff_s_inner(j,i), nzt_diff
	550	tend(k,j,i) = tend(k,j,i) - &
	551	kh(k,j,i) * g / pt_reference * &
	552	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
	553	dd2zu(k)
	554	ENDDO
[97]	555
[940]	556	IF ( use_surface_fluxes ) THEN
	557	k = nzb_diff_s_inner(j,i)-1
	558	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
	559	shf(j,i)
	560	ENDIF
[97]	561
[940]	562	IF ( use_top_fluxes ) THEN
	563	k = nzt
	564	tend(k,j,i) = tend(k,j,i) + g / pt_reference * &
	565	tswst(j,i)
	566	ENDIF
	567	ENDDO
[57]	568
[940]	569	ENDIF
[57]	570
[940]	571	ELSE
[1]	572
[940]	573	IF ( ocean ) THEN
[97]	574	!
[940]	575	!-- So far in the ocean no special treatment of density flux
	576	!-- in the bottom and top surface layer
	577	DO j = nys, nyn
	578	DO k = nzb_s_inner(j,i)+1, nzt
	579	tend(k,j,i) = tend(k,j,i) + &
[2031]	580	kh(k,j,i) * g / rho_ocean(k,j,i) * &
	581	( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * &
[940]	582	dd2zu(k)
	583	ENDDO
[97]	584	ENDDO
	585
[940]	586	ELSE
[97]	587
[940]	588	DO j = nys, nyn
	589	DO k = nzb_diff_s_inner(j,i), nzt_diff
	590	tend(k,j,i) = tend(k,j,i) - &
	591	kh(k,j,i) * g / pt(k,j,i) * &
	592	( pt(k+1,j,i) - pt(k-1,j,i) ) * &
	593	dd2zu(k)
	594	ENDDO
	595
	596	IF ( use_surface_fluxes ) THEN
	597	k = nzb_diff_s_inner(j,i)-1
	598	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
	599	shf(j,i)
	600	ENDIF
	601
	602	IF ( use_top_fluxes ) THEN
	603	k = nzt
	604	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * &
	605	tswst(j,i)
	606	ENDIF
[97]	607	ENDDO
	608
[940]	609	ENDIF
[97]	610
	611	ENDIF
[1]	612
[940]	613	ELSE
[57]	614
[940]	615	DO j = nys, nyn
[1]	616
[940]	617	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	618
[1007]	619	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	620	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	621	k2 = 0.61_wp * pt(k,j,i)
[1007]	622	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	623	g / vpt(k,j,i) * &
	624	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	625	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	626	) * dd2zu(k)
	627	ELSE IF ( cloud_physics ) THEN
[1342]	628	IF ( ql(k,j,i) == 0.0_wp ) THEN
	629	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	630	k2 = 0.61_wp * pt(k,j,i)
[940]	631	ELSE
	632	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	633	temp = theta * t_d_pt(k)
[1342]	634	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	635	( q(k,j,i) - ql(k,j,i) ) * &
	636	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	637	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	638	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	639	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	640	ENDIF
[1007]	641	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * &
	642	g / vpt(k,j,i) * &
[940]	643	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	644	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	645	) * dd2zu(k)
[1007]	646	ELSE IF ( cloud_droplets ) THEN
[1342]	647	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	648	k2 = 0.61_wp * pt(k,j,i)
[1007]	649	tend(k,j,i) = tend(k,j,i) - &
	650	kh(k,j,i) * g / vpt(k,j,i) * &
	651	( k1 * ( pt(k+1,j,i)- pt(k-1,j,i) ) + &
	652	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
	653	pt(k,j,i) * ( ql(k+1,j,i) - &
	654	ql(k-1,j,i) ) ) * dd2zu(k)
	655	ENDIF
	656
[940]	657	ENDDO
	658
[1]	659	ENDDO
	660
[940]	661	IF ( use_surface_fluxes ) THEN
[1]	662
[940]	663	DO j = nys, nyn
[1]	664
[940]	665	k = nzb_diff_s_inner(j,i)-1
[1]	666
[1007]	667	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	668	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	669	k2 = 0.61_wp * pt(k,j,i)
[1007]	670	ELSE IF ( cloud_physics ) THEN
[1342]	671	IF ( ql(k,j,i) == 0.0_wp ) THEN
	672	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	673	k2 = 0.61_wp * pt(k,j,i)
[940]	674	ELSE
	675	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	676	temp = theta * t_d_pt(k)
[1342]	677	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
[1353]	678	( q(k,j,i) - ql(k,j,i) ) * &
	679	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	680	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	681	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	682	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	683	ENDIF
[1007]	684	ELSE IF ( cloud_droplets ) THEN
[1342]	685	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	686	k2 = 0.61_wp * pt(k,j,i)
[1]	687	ENDIF
	688
[940]	689	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	690	( k1* shf(j,i) + k2 * qsws(j,i) )
	691	ENDDO
[1]	692
[940]	693	ENDIF
[1]	694
[940]	695	IF ( use_top_fluxes ) THEN
[19]	696
[940]	697	DO j = nys, nyn
[19]	698
[940]	699	k = nzt
[19]	700
[1007]	701	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	702	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	703	k2 = 0.61_wp * pt(k,j,i)
[1007]	704	ELSE IF ( cloud_physics ) THEN
[1342]	705	IF ( ql(k,j,i) == 0.0_wp ) THEN
	706	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	707	k2 = 0.61_wp * pt(k,j,i)
[940]	708	ELSE
	709	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	710	temp = theta * t_d_pt(k)
[1353]	711	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	712	( q(k,j,i) - ql(k,j,i) ) * &
	713	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	714	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	715	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	716	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	717	ENDIF
[1007]	718	ELSE IF ( cloud_droplets ) THEN
[1342]	719	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	720	k2 = 0.61_wp * pt(k,j,i)
[19]	721	ENDIF
	722
[940]	723	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	724	( k1* tswst(j,i) + k2 * qswst(j,i) )
	725	ENDDO
[19]	726
[940]	727	ENDIF
	728
[19]	729	ENDIF
	730
[1]	731	ENDIF
	732
	733	ENDDO
	734
	735	END SUBROUTINE production_e
	736
	737
	738	!------------------------------------------------------------------------------!
[1682]	739	! Description:
	740	! ------------
	741	!> Call for grid point i,j
[1]	742	!------------------------------------------------------------------------------!
	743	SUBROUTINE production_e_ij( i, j )
	744
[1320]	745	USE arrays_3d, &
[2031]	746	ONLY: ddzw, dd2zu, kh, km, pt, q, ql, qsws, qswst, rho_ocean, shf, &
[1320]	747	tend, tswst, u, v, vpt, w
[449]	748
[1320]	749	USE cloud_parameters, &
	750	ONLY: l_d_cp, l_d_r, pt_d_t, t_d_pt
	751
	752	USE control_parameters, &
[1691]	753	ONLY: cloud_droplets, cloud_physics, constant_flux_layer, g, &
	754	humidity, kappa, neutral, ocean, pt_reference, &
	755	rho_reference, use_single_reference_value, &
	756	use_surface_fluxes, use_top_fluxes
[1320]	757
	758	USE grid_variables, &
	759	ONLY: ddx, dx, ddy, dy, wall_e_x, wall_e_y
	760
	761	USE indices, &
	762	ONLY: nxl, nxr, nys, nyn, nzb, nzb_diff_s_inner, &
	763	nzb_diff_s_outer, nzb_s_inner, nzt, nzt_diff
	764
[1]	765	IMPLICIT NONE
	766
[1682]	767	INTEGER(iwp) :: i !<
	768	INTEGER(iwp) :: j !<
	769	INTEGER(iwp) :: k !<
[1]	770
[1682]	771	REAL(wp) :: def !<
	772	REAL(wp) :: dudx !<
	773	REAL(wp) :: dudy !<
	774	REAL(wp) :: dudz !<
	775	REAL(wp) :: dvdx !<
	776	REAL(wp) :: dvdy !<
	777	REAL(wp) :: dvdz !<
	778	REAL(wp) :: dwdx !<
	779	REAL(wp) :: dwdy !<
	780	REAL(wp) :: dwdz !<
	781	REAL(wp) :: k1 !<
	782	REAL(wp) :: k2 !<
	783	REAL(wp) :: km_neutral !<
	784	REAL(wp) :: theta !<
	785	REAL(wp) :: temp !<
[1]	786
[1682]	787	REAL(wp), DIMENSION(nzb:nzt+1) :: usvs !<
	788	REAL(wp), DIMENSION(nzb:nzt+1) :: vsus !<
	789	REAL(wp), DIMENSION(nzb:nzt+1) :: wsus !<
	790	REAL(wp), DIMENSION(nzb:nzt+1) :: wsvs !<
[53]	791
[1]	792	!
	793	!-- Calculate TKE production by shear
[19]	794	DO k = nzb_diff_s_outer(j,i), nzt
[1]	795
[1342]	796	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	797	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	798	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	799	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	800	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[1]	801
[1342]	802	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	803	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	804	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	805	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	806	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[1]	807
[1342]	808	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	809	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	810	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	811	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	812	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	813
[1342]	814	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) &
	815	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
	816	+ 2.0_wp * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	817
[1342]	818	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	819
	820	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
[1007]	821
[1]	822	ENDDO
	823
[1691]	824	IF ( constant_flux_layer ) THEN
[1]	825
[1342]	826	IF ( ( wall_e_x(j,i) /= 0.0_wp ) .OR. ( wall_e_y(j,i) /= 0.0_wp ) ) THEN
[55]	827
[1]	828	!
[55]	829	!-- Position beneath wall
	830	!-- (2) - Will allways be executed.
	831	!-- 'bottom and wall: use u_0,v_0 and wall functions'
[1]	832	k = nzb_diff_s_inner(j,i)-1
	833
	834	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[1342]	835	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	836	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
	837	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	838	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	839	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
[53]	840	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	841
[1342]	842	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
[1007]	843	!
[208]	844	!-- Inconsistency removed: as the thermal stratification
	845	!-- is not taken into account for the evaluation of the
	846	!-- wall fluxes at vertical walls, the eddy viscosity km
	847	!-- must not be used for the evaluation of the velocity
	848	!-- gradients dudy and dwdy
	849	!-- Note: The validity of the new method has not yet
	850	!-- been shown, as so far no suitable data for a
	851	!-- validation has been available
[53]	852	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	853	usvs, 1.0_wp, 0.0_wp, 0.0_wp, 0.0_wp )
[53]	854	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	855	wsvs, 0.0_wp, 0.0_wp, 1.0_wp, 0.0_wp )
[1342]	856	km_neutral = kappa * ( usvs(k)2 + wsvs(k)2 )*0.25_wp &
	857	0.5_wp * dy
	858	IF ( km_neutral > 0.0_wp ) THEN
[364]	859	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	860	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	861	ELSE
[1342]	862	dudy = 0.0_wp
	863	dwdy = 0.0_wp
[364]	864	ENDIF
[1]	865	ELSE
[1342]	866	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	867	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	868	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	869	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	870	ENDIF
	871
[1342]	872	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
[1007]	873	!
[208]	874	!-- Inconsistency removed: as the thermal stratification
	875	!-- is not taken into account for the evaluation of the
	876	!-- wall fluxes at vertical walls, the eddy viscosity km
	877	!-- must not be used for the evaluation of the velocity
	878	!-- gradients dvdx and dwdx
	879	!-- Note: The validity of the new method has not yet
	880	!-- been shown, as so far no suitable data for a
	881	!-- validation has been available
[53]	882	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	883	vsus, 0.0_wp, 1.0_wp, 0.0_wp, 0.0_wp )
[53]	884	CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
[1320]	885	wsus, 0.0_wp, 0.0_wp, 0.0_wp, 1.0_wp )
[1342]	886	km_neutral = kappa * ( vsus(k)2 + wsus(k)2 )*0.25_wp &
	887	0.5_wp * dx
	888	IF ( km_neutral > 0.0_wp ) THEN
[364]	889	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	890	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	891	ELSE
[1342]	892	dvdx = 0.0_wp
	893	dwdx = 0.0_wp
[364]	894	ENDIF
[1]	895	ELSE
[1342]	896	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	897	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	898	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	899	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
[1]	900	ENDIF
	901
[1342]	902	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	903	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	904	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	905
[1342]	906	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	907
	908	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	909
	910	!
[55]	911	!-- (3) - will be executed only, if there is at least one level
	912	!-- between (2) and (4), i.e. the topography must have a
	913	!-- minimum height of 2 dz. Wall fluxes for this case have
	914	!-- already been calculated for (2).
	915	!-- 'wall only: use wall functions'
[1]	916	DO k = nzb_diff_s_inner(j,i), nzb_diff_s_outer(j,i)-2
	917
	918	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
[1342]	919	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	920	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
	921	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	922	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	923	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[53]	924	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	925
[1342]	926	IF ( wall_e_y(j,i) /= 0.0_wp ) THEN
[1007]	927	!
[208]	928	!-- Inconsistency removed: as the thermal stratification
	929	!-- is not taken into account for the evaluation of the
	930	!-- wall fluxes at vertical walls, the eddy viscosity km
	931	!-- must not be used for the evaluation of the velocity
	932	!-- gradients dudy and dwdy
	933	!-- Note: The validity of the new method has not yet
	934	!-- been shown, as so far no suitable data for a
	935	!-- validation has been available
	936	km_neutral = kappa * ( usvs(k)**2 + &
[1342]	937	wsvs(k)2 )0.25_wp * 0.5_wp * dy
	938	IF ( km_neutral > 0.0_wp ) THEN
[364]	939	dudy = - wall_e_y(j,i) * usvs(k) / km_neutral
	940	dwdy = - wall_e_y(j,i) * wsvs(k) / km_neutral
	941	ELSE
[1342]	942	dudy = 0.0_wp
	943	dwdy = 0.0_wp
[364]	944	ENDIF
[1]	945	ELSE
[1342]	946	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	947	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	948	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	949	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
[1]	950	ENDIF
	951
[1342]	952	IF ( wall_e_x(j,i) /= 0.0_wp ) THEN
[1007]	953	!
[208]	954	!-- Inconsistency removed: as the thermal stratification
	955	!-- is not taken into account for the evaluation of the
	956	!-- wall fluxes at vertical walls, the eddy viscosity km
	957	!-- must not be used for the evaluation of the velocity
	958	!-- gradients dvdx and dwdx
	959	!-- Note: The validity of the new method has not yet
	960	!-- been shown, as so far no suitable data for a
	961	!-- validation has been available
	962	km_neutral = kappa * ( vsus(k)**2 + &
[1342]	963	wsus(k)2 )0.25_wp * 0.5_wp * dx
	964	IF ( km_neutral > 0.0_wp ) THEN
[364]	965	dvdx = - wall_e_x(j,i) * vsus(k) / km_neutral
	966	dwdx = - wall_e_x(j,i) * wsus(k) / km_neutral
	967	ELSE
[1342]	968	dvdx = 0.0_wp
	969	dwdx = 0.0_wp
[364]	970	ENDIF
[1]	971	ELSE
[1342]	972	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	973	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	974	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	975	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
[1]	976	ENDIF
	977
[1342]	978	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	979	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	980	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	981
[1342]	982	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	983
	984	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	985
	986	ENDDO
	987
	988	!
[55]	989	!-- (4) - will allways be executed.
	990	!-- 'special case: free atmosphere' (as for case (0))
[1]	991	k = nzb_diff_s_outer(j,i)-1
	992
[1342]	993	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	994	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	995	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	996	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	997	u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[1]	998
[1342]	999	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1000	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1001	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1002	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1003	v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[1]	1004
[1342]	1005	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1006	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1007	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1008	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1009	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	1010
[1353]	1011	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[1]	1012	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1353]	1013	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	1014
[1342]	1015	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	1016
	1017	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1018
	1019	ELSE
	1020
	1021	!
[55]	1022	!-- Position without adjacent wall
	1023	!-- (1) - will allways be executed.
	1024	!-- 'bottom only: use u_0,v_0'
[1]	1025	k = nzb_diff_s_inner(j,i)-1
	1026
[1342]	1027	dudx = ( u(k,j,i+1) - u(k,j,i) ) * ddx
	1028	dudy = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - &
	1029	u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
	1030	dudz = 0.5_wp * ( u(k+1,j,i) + u(k+1,j,i+1) - &
	1031	u_0(j,i) - u_0(j,i+1) ) * dd2zu(k)
[1]	1032
[1342]	1033	dvdx = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - &
	1034	v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
	1035	dvdy = ( v(k,j+1,i) - v(k,j,i) ) * ddy
	1036	dvdz = 0.5_wp * ( v(k+1,j,i) + v(k+1,j+1,i) - &
	1037	v_0(j,i) - v_0(j+1,i) ) * dd2zu(k)
[1]	1038
[1342]	1039	dwdx = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - &
	1040	w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
	1041	dwdy = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - &
	1042	w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
	1043	dwdz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[1]	1044
[1342]	1045	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) &
[1]	1046	+ dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz2 + dvdz2 &
[1342]	1047	+ 2.0_wp * ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[1]	1048
[1342]	1049	IF ( def < 0.0_wp ) def = 0.0_wp
[1]	1050
	1051	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1052
	1053	ENDIF
	1054
[37]	1055	ELSEIF ( use_surface_fluxes ) THEN
	1056
	1057	k = nzb_diff_s_outer(j,i)-1
	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(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
[37]	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(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
[37]	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)
[37]	1076
[1342]	1077	def = 2.0_wp * ( dudx2 + dvdy2 + dwdz**2 ) + &
[37]	1078	dudy2 + dvdx2 + dwdx2 + dwdy2 + dudz**2 + &
[1342]	1079	dvdz*2 + 2.0_wp ( dvdxdudy + dwdxdudz + dwdy*dvdz )
[37]	1080
[1342]	1081	IF ( def < 0.0_wp ) def = 0.0_wp
[37]	1082
	1083	tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def
	1084
[1]	1085	ENDIF
	1086
	1087	!
[940]	1088	!-- If required, calculate TKE production by buoyancy
	1089	IF ( .NOT. neutral ) THEN
[1]	1090
[940]	1091	IF ( .NOT. humidity ) THEN
[19]	1092
[1179]	1093	IF ( use_single_reference_value ) THEN
[940]	1094
	1095	IF ( ocean ) THEN
[97]	1096	!
[940]	1097	!-- So far in the ocean no special treatment of density flux in
	1098	!-- the bottom and top surface layer
	1099	DO k = nzb_s_inner(j,i)+1, nzt
	1100	tend(k,j,i) = tend(k,j,i) + &
	1101	kh(k,j,i) * g / rho_reference * &
[2031]	1102	( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * dd2zu(k)
[940]	1103	ENDDO
[97]	1104
[940]	1105	ELSE
[97]	1106
[940]	1107	DO k = nzb_diff_s_inner(j,i), nzt_diff
	1108	tend(k,j,i) = tend(k,j,i) - &
	1109	kh(k,j,i) * g / pt_reference * &
	1110	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	1111	ENDDO
[1]	1112
[940]	1113	IF ( use_surface_fluxes ) THEN
	1114	k = nzb_diff_s_inner(j,i)-1
	1115	tend(k,j,i) = tend(k,j,i) + g / pt_reference * shf(j,i)
	1116	ENDIF
[19]	1117
[940]	1118	IF ( use_top_fluxes ) THEN
	1119	k = nzt
	1120	tend(k,j,i) = tend(k,j,i) + g / pt_reference * tswst(j,i)
	1121	ENDIF
	1122
[97]	1123	ENDIF
	1124
[940]	1125	ELSE
[57]	1126
[940]	1127	IF ( ocean ) THEN
[97]	1128	!
[940]	1129	!-- So far in the ocean no special treatment of density flux in
	1130	!-- the bottom and top surface layer
	1131	DO k = nzb_s_inner(j,i)+1, nzt
	1132	tend(k,j,i) = tend(k,j,i) + &
[2031]	1133	kh(k,j,i) * g / rho_ocean(k,j,i) * &
	1134	( rho_ocean(k+1,j,i) - rho_ocean(k-1,j,i) ) * dd2zu(k)
[940]	1135	ENDDO
[97]	1136
[940]	1137	ELSE
[97]	1138
[940]	1139	DO k = nzb_diff_s_inner(j,i), nzt_diff
	1140	tend(k,j,i) = tend(k,j,i) - &
	1141	kh(k,j,i) * g / pt(k,j,i) * &
	1142	( pt(k+1,j,i) - pt(k-1,j,i) ) * dd2zu(k)
	1143	ENDDO
[57]	1144
[940]	1145	IF ( use_surface_fluxes ) THEN
	1146	k = nzb_diff_s_inner(j,i)-1
	1147	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * shf(j,i)
	1148	ENDIF
[57]	1149
[940]	1150	IF ( use_top_fluxes ) THEN
	1151	k = nzt
	1152	tend(k,j,i) = tend(k,j,i) + g / pt(k,j,i) * tswst(j,i)
	1153	ENDIF
	1154
[97]	1155	ENDIF
	1156
[57]	1157	ENDIF
	1158
[940]	1159	ELSE
[57]	1160
[940]	1161	DO k = nzb_diff_s_inner(j,i), nzt_diff
[1]	1162
[1007]	1163	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1164	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1165	k2 = 0.61_wp * pt(k,j,i)
[1007]	1166	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1167	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1168	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1169	) * dd2zu(k)
	1170	ELSE IF ( cloud_physics ) THEN
[1342]	1171	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1172	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1173	k2 = 0.61_wp * pt(k,j,i)
[940]	1174	ELSE
	1175	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1176	temp = theta * t_d_pt(k)
[1342]	1177	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1178	( q(k,j,i) - ql(k,j,i) ) * &
	1179	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1180	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	1181	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	1182	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	1183	ENDIF
[1007]	1184	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
[940]	1185	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1186	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) &
	1187	) * dd2zu(k)
[1007]	1188	ELSE IF ( cloud_droplets ) THEN
[1342]	1189	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1190	k2 = 0.61_wp * pt(k,j,i)
[1007]	1191	tend(k,j,i) = tend(k,j,i) - kh(k,j,i) * g / vpt(k,j,i) * &
	1192	( k1 * ( pt(k+1,j,i)-pt(k-1,j,i) ) + &
	1193	k2 * ( q(k+1,j,i) - q(k-1,j,i) ) - &
	1194	pt(k,j,i) * ( ql(k+1,j,i) - &
	1195	ql(k-1,j,i) ) ) * dd2zu(k)
	1196	ENDIF
[940]	1197	ENDDO
[19]	1198
[940]	1199	IF ( use_surface_fluxes ) THEN
	1200	k = nzb_diff_s_inner(j,i)-1
[1]	1201
[1007]	1202	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1203	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1204	k2 = 0.61_wp * pt(k,j,i)
[1007]	1205	ELSE IF ( cloud_physics ) THEN
[1342]	1206	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1207	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1208	k2 = 0.61_wp * pt(k,j,i)
[940]	1209	ELSE
	1210	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1211	temp = theta * t_d_pt(k)
[1342]	1212	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1213	( q(k,j,i) - ql(k,j,i) ) * &
	1214	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1215	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	1216	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	1217	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	1218	ENDIF
[1007]	1219	ELSE IF ( cloud_droplets ) THEN
[1342]	1220	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1221	k2 = 0.61_wp * pt(k,j,i)
[1]	1222	ENDIF
[940]	1223
	1224	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1225	( k1* shf(j,i) + k2 * qsws(j,i) )
[1]	1226	ENDIF
	1227
[940]	1228	IF ( use_top_fluxes ) THEN
	1229	k = nzt
[1]	1230
[1007]	1231	IF ( .NOT. cloud_physics .AND. .NOT. cloud_droplets ) THEN
[1342]	1232	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1233	k2 = 0.61_wp * pt(k,j,i)
[1007]	1234	ELSE IF ( cloud_physics ) THEN
[1342]	1235	IF ( ql(k,j,i) == 0.0_wp ) THEN
	1236	k1 = 1.0_wp + 0.61_wp * q(k,j,i)
	1237	k2 = 0.61_wp * pt(k,j,i)
[940]	1238	ELSE
	1239	theta = pt(k,j,i) + pt_d_t(k) * l_d_cp * ql(k,j,i)
	1240	temp = theta * t_d_pt(k)
[1342]	1241	k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * &
	1242	( q(k,j,i) - ql(k,j,i) ) * &
	1243	( 1.0_wp + 0.622_wp * l_d_r / temp ) ) / &
	1244	( 1.0_wp + 0.622_wp * l_d_r * l_d_cp * &
[940]	1245	( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
[1342]	1246	k2 = theta * ( l_d_cp / temp * k1 - 1.0_wp )
[940]	1247	ENDIF
[1007]	1248	ELSE IF ( cloud_droplets ) THEN
[1342]	1249	k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
	1250	k2 = 0.61_wp * pt(k,j,i)
[19]	1251	ENDIF
[940]	1252
	1253	tend(k,j,i) = tend(k,j,i) + g / vpt(k,j,i) * &
	1254	( k1* tswst(j,i) + k2 * qswst(j,i) )
[19]	1255	ENDIF
	1256
	1257	ENDIF
	1258
[1]	1259	ENDIF
	1260
	1261	END SUBROUTINE production_e_ij
	1262
	1263
[1682]	1264	!------------------------------------------------------------------------------!
	1265	! Description:
	1266	! ------------
	1267	!> @todo Missing subroutine description.
	1268	!------------------------------------------------------------------------------!
[1]	1269	SUBROUTINE production_e_init
	1270
[1320]	1271	USE arrays_3d, &
	1272	ONLY: kh, km, u, us, usws, v, vsws, zu
[1]	1273
[1320]	1274	USE control_parameters, &
[1691]	1275	ONLY: constant_flux_layer, kappa
[1320]	1276
	1277	USE indices, &
	1278	ONLY: nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb_u_inner, &
	1279	nzb_v_inner
	1280
[1]	1281	IMPLICIT NONE
	1282
[1682]	1283	INTEGER(iwp) :: i !<
	1284	INTEGER(iwp) :: j !<
	1285	INTEGER(iwp) :: ku !<
	1286	INTEGER(iwp) :: kv !<
[1]	1287
[1691]	1288	IF ( constant_flux_layer ) THEN
[1]	1289
	1290	IF ( first_call ) THEN
[759]	1291	ALLOCATE( u_0(nysg:nyng,nxlg:nxrg), v_0(nysg:nyng,nxlg:nxrg) )
[1342]	1292	u_0 = 0.0_wp ! just to avoid access of uninitialized memory
	1293	v_0 = 0.0_wp ! within exchange_horiz_2d
[1]	1294	first_call = .FALSE.
	1295	ENDIF
	1296
	1297	!
	1298	!-- Calculate a virtual velocity at the surface in a way that the
	1299	!-- vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
	1300	!-- Prandtl law (-w'u'/km). This gradient is used in the TKE shear
	1301	!-- production term at k=1 (see production_e_ij).
	1302	!-- The velocity gradient has to be limited in case of too small km
	1303	!-- (otherwise the timestep may be significantly reduced by large
	1304	!-- surface winds).
[106]	1305	!-- Upper bounds are nxr+1 and nyn+1 because otherwise these values are
	1306	!-- not available in case of non-cyclic boundary conditions.
[1]	1307	!-- WARNING: the exact analytical solution would require the determination
	1308	!-- of the eddy diffusivity by km = u* * kappa * zp / phi_m.
	1309	!$OMP PARALLEL DO PRIVATE( ku, kv )
[106]	1310	DO i = nxl, nxr+1
	1311	DO j = nys, nyn+1
[1]	1312
	1313	ku = nzb_u_inner(j,i)+1
	1314	kv = nzb_v_inner(j,i)+1
	1315
	1316	u_0(j,i) = u(ku+1,j,i) + usws(j,i) * ( zu(ku+1) - zu(ku-1) ) / &
[1342]	1317	( 0.5_wp * ( km(ku,j,i) + km(ku,j,i-1) ) + &
	1318	1.0E-20_wp )
[1]	1319	! ( us(j,i) * kappa * zu(1) )
	1320	v_0(j,i) = v(kv+1,j,i) + vsws(j,i) * ( zu(kv+1) - zu(kv-1) ) / &
[1342]	1321	( 0.5_wp * ( km(kv,j,i) + km(kv,j-1,i) ) + &
	1322	1.0E-20_wp )
[1]	1323	! ( us(j,i) * kappa * zu(1) )
	1324
	1325	IF ( ABS( u(ku+1,j,i) - u_0(j,i) ) > &
	1326	ABS( u(ku+1,j,i) - u(ku-1,j,i) ) ) u_0(j,i) = u(ku-1,j,i)
	1327	IF ( ABS( v(kv+1,j,i) - v_0(j,i) ) > &
	1328	ABS( v(kv+1,j,i) - v(kv-1,j,i) ) ) v_0(j,i) = v(kv-1,j,i)
	1329
	1330	ENDDO
	1331	ENDDO
	1332
	1333	CALL exchange_horiz_2d( u_0 )
	1334	CALL exchange_horiz_2d( v_0 )
	1335
	1336	ENDIF
	1337
	1338	END SUBROUTINE production_e_init
	1339
	1340	END MODULE production_e_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |