Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |