Home

Context Navigation

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

Last change on this file since 1016 was 1015, checked in by raasch, 12 years ago
Starting with changes required for GPU optimization. OpenACC statements for using NVIDIA GPUs added. Adjustment of mixing length to the Prandtl mixing length at first grid point above ground removed. mask array is set zero for ghost boundaries
Property svn:keywords set to `Id`
File size: 82.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |