Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |