Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |