Home

Context Navigation

advec_ws.f90 @ 1221

Last change on this file since 1221 was 1221, checked in by raasch, 12 years ago

New:

openACC porting of reduction operations
additional 3D-flag arrays for replacing the 2D-index arrays nzb_s_inner and nzb_diff_s_inner
(flow_statistics, init_grid, init_3d_model, modules, palm, pres, time_integration)

Changed:

for PGI/openACC performance reasons set default compile options for openACC to "-ta=nocache",
and set environment variable PGI_ACC_SYNCHRONOUS=1
(MAKE.inc.pgi.openacc, palm_simple_run)

wall_flags_0 changed to 32bit INTEGER, additional array wall_flags_00 introduced to hold
bits 32-63
(advec_ws, init_grid, modules, palm)

Errors:

dummy argument tri in 1d-routines replaced by tri_for_1d because of name
conflict with arry tri in module arrays_3d
(tridia_solver)

Property svn:executable set to *
Property svn:keywords set to Id

File size: 272.0 KB

Rev	Line
[713]	1	MODULE advec_ws
	2
[1036]	3	!--------------------------------------------------------------------------------!
	4	! This file is part of PALM.
	5	!
	6	! PALM is free software: you can redistribute it and/or modify it under the terms
	7	! of the GNU General Public License as published by the Free Software Foundation,
	8	! either version 3 of the License, or (at your option) any later version.
	9	!
	10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	13	!
	14	! You should have received a copy of the GNU General Public License along with
	15	! PALM. If not, see <http://www.gnu.org/licenses/>.
	16	!
	17	! Copyright 1997-2012 Leibniz University Hannover
	18	!--------------------------------------------------------------------------------!
	19	!
[713]	20	! Current revisions:
[736]	21	! ------------------
[1221]	22	! wall_flags_00 introduced, which holds bits 32-...
[863]	23	!
	24	! Former revisions:
	25	! -----------------
	26	! $Id: advec_ws.f90 1221 2013-09-10 08:59:13Z raasch $
[889]	27	!
[1132]	28	! 1128 2013-04-12 06:19:32Z raasch
	29	! loop index bounds in accelerator version replaced by i_left, i_right, j_south,
	30	! j_north
	31	!
[1116]	32	! 1115 2013-03-26 18:16:16Z hoffmann
	33	! calculation of qr and nr is restricted to precipitation
	34	!
[1054]	35	! 1053 2012-11-13 17:11:03Z hoffmann
	36	! necessary expansions according to the two new prognostic equations (nr, qr)
	37	! of the two-moment cloud physics scheme:
	38	! +flux_l_, flux_s_, diss_l_, diss_s_, sums_ws*s_ws_l
	39	!
[1037]	40	! 1036 2012-10-22 13:43:42Z raasch
	41	! code put under GPL (PALM 3.9)
	42	!
[1028]	43	! 1027 2012-10-15 17:18:39Z suehring
	44	! Bugfix in calculation indices k_mm, k_pp in accelerator version
	45	!
[1020]	46	! 1019 2012-09-28 06:46:45Z raasch
	47	! small change in comment lines
	48	!
[1017]	49	! 1015 2012-09-27 09:23:24Z raasch
	50	! accelerator versions (*_acc) added
	51	!
[1011]	52	! 1010 2012-09-20 07:59:54Z raasch
	53	! cpp switch __nopointer added for pointer free version
	54	!
[889]	55	! 888 2012-04-20 15:03:46Z suehring
	56	! Number of IBITS() calls with identical arguments is reduced.
	57	!
[863]	58	! 862 2012-03-26 14:21:38Z suehring
[861]	59	! ws-scheme also work with topography in combination with vector version.
	60	! ws-scheme also work with outflow boundaries in combination with
	61	! vector version.
	62	! Degradation of the applied order of scheme is now steered by multiplying with
	63	! Integer wall_flags_0. 2nd order scheme, WS3 and WS5 are calculated on each
	64	! grid point and mulitplied with the appropriate flag.
	65	! 2nd order numerical dissipation term changed. Now the appropriate 2nd order
	66	! term derived according to the 4th and 6th order terms is applied. It turns
	67	! out that diss_2nd does not provide sufficient dissipation near walls.
	68	! Therefore, the function diss_2nd is removed.
	69	! Near walls a divergence correction is necessary to overcome numerical
	70	! instabilities due to too less divergence reduction of the velocity field.
	71	! boundary_flags and logicals steering the degradation are removed.
	72	! Empty SUBROUTINE local_diss removed.
[862]	73	! Further formatting adjustments.
	74	!
[802]	75	! 801 2012-01-10 17:30:36Z suehring
[801]	76	! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
	77	! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
	78	! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
	79	! dimension.
[737]	80	!
[744]	81	! 743 2011-08-18 16:10:16Z suehring
	82	! Evaluation of turbulent fluxes with WS-scheme only for the whole model
	83	! domain. Therefor dimension of arrays needed for statistical evaluation
	84	! decreased.
	85	!
[737]	86	! 736 2011-08-17 14:13:26Z suehring
[736]	87	! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
	88	! index where fluxes on left side have to be calculated explicitly is
	89	! different on each thread. Furthermore the swapping of fluxes is now
	90	! thread-safe by adding an additional dimension.
[714]	91	!
	92	! 713 2011-03-30 14:21:21Z suehring
[713]	93	! File reformatted.
	94	! Bugfix in vertical advection of w concerning the optimized version for
	95	! vector architecture.
	96	! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
	97	! broken numbers.
	98	!
	99	! 709 2011-03-30 09:31:40Z raasch
	100	! formatting adjustments
	101	!
	102	! 705 2011-03-25 11:21:43 Z suehring $
	103	! Bugfix in declaration of logicals concerning outflow boundaries.
	104	!
	105	! 411 2009-12-11 12:31:43 Z suehring
	106	! Allocation of weight_substep moved to init_3d_model.
	107	! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters.
	108	! Setting bc for the horizontal velocity variances added (moved from
	109	! flow_statistics).
	110	!
	111	! Initial revision
	112	!
	113	! 411 2009-12-11 12:31:43 Z suehring
	114	!
	115	! Description:
	116	! ------------
	117	! Advection scheme for scalars and momentum using the flux formulation of
	118	! Wicker and Skamarock 5th order. Additionally the module contains of a
	119	! routine using for initialisation and steering of the statical evaluation.
	120	! The computation of turbulent fluxes takes place inside the advection
[861]	121	! routines.
	122	! Near non-cyclic boundaries the order of the applied advection scheme is
	123	! degraded.
	124	! A divergence correction is applied. It is necessary for topography, since
	125	! the divergence is not sufficiently reduced, resulting in erroneous fluxes and
	126	! partly numerical instabilities.
[713]	127	!-----------------------------------------------------------------------------!
	128
	129	PRIVATE
[1015]	130	PUBLIC advec_s_ws, advec_s_ws_acc, advec_u_ws, advec_u_ws_acc, &
	131	advec_v_ws, advec_v_ws_acc, advec_w_ws, advec_w_ws_acc, &
[861]	132	ws_init, ws_statistics
[713]	133
	134	INTERFACE ws_init
	135	MODULE PROCEDURE ws_init
	136	END INTERFACE ws_init
	137
	138	INTERFACE ws_statistics
	139	MODULE PROCEDURE ws_statistics
	140	END INTERFACE ws_statistics
	141
	142	INTERFACE advec_s_ws
	143	MODULE PROCEDURE advec_s_ws
	144	MODULE PROCEDURE advec_s_ws_ij
	145	END INTERFACE advec_s_ws
	146
	147	INTERFACE advec_u_ws
	148	MODULE PROCEDURE advec_u_ws
	149	MODULE PROCEDURE advec_u_ws_ij
	150	END INTERFACE advec_u_ws
	151
[1015]	152	INTERFACE advec_u_ws_acc
	153	MODULE PROCEDURE advec_u_ws_acc
	154	END INTERFACE advec_u_ws_acc
	155
[713]	156	INTERFACE advec_v_ws
	157	MODULE PROCEDURE advec_v_ws
	158	MODULE PROCEDURE advec_v_ws_ij
	159	END INTERFACE advec_v_ws
	160
[1015]	161	INTERFACE advec_v_ws_acc
	162	MODULE PROCEDURE advec_v_ws_acc
	163	END INTERFACE advec_v_ws_acc
	164
[713]	165	INTERFACE advec_w_ws
	166	MODULE PROCEDURE advec_w_ws
	167	MODULE PROCEDURE advec_w_ws_ij
	168	END INTERFACE advec_w_ws
	169
[1015]	170	INTERFACE advec_w_ws_acc
	171	MODULE PROCEDURE advec_w_ws_acc
	172	END INTERFACE advec_w_ws_acc
	173
[713]	174	CONTAINS
	175
	176
	177	!------------------------------------------------------------------------------!
	178	! Initialization of WS-scheme
	179	!------------------------------------------------------------------------------!
	180	SUBROUTINE ws_init
	181
	182	USE arrays_3d
	183	USE constants
	184	USE control_parameters
	185	USE indices
[736]	186	USE pegrid
[713]	187	USE statistics
	188
	189	!
	190	!-- Set the appropriate factors for scalar and momentum advection.
	191	adv_sca_5 = 1./60.
	192	adv_sca_3 = 1./12.
[861]	193	adv_sca_1 = 1./2.
[713]	194	adv_mom_5 = 1./120.
	195	adv_mom_3 = 1./24.
[861]	196	adv_mom_1 = 1./4.
[713]	197	!
	198	!-- Arrays needed for statical evaluation of fluxes.
	199	IF ( ws_scheme_mom ) THEN
	200
[801]	201	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
	202	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
	203	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
	204	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
	205	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
[713]	206
	207	sums_wsus_ws_l = 0.0
	208	sums_wsvs_ws_l = 0.0
	209	sums_us2_ws_l = 0.0
	210	sums_vs2_ws_l = 0.0
	211	sums_ws2_ws_l = 0.0
	212
	213	ENDIF
	214
	215	IF ( ws_scheme_sca ) THEN
	216
[801]	217	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
[713]	218	sums_wspts_ws_l = 0.0
	219
[1115]	220	IF ( humidity .OR. passive_scalar ) THEN
[801]	221	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
[713]	222	sums_wsqs_ws_l = 0.0
	223	ENDIF
	224
[1115]	225	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	226	precipitation ) THEN
[1053]	227	ALLOCATE( sums_wsqrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
	228	ALLOCATE( sums_wsnrs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
	229	sums_wsqrs_ws_l = 0.0
	230	sums_wsnrs_ws_l = 0.0
	231	ENDIF
	232
[713]	233	IF ( ocean ) THEN
[801]	234	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
[713]	235	sums_wssas_ws_l = 0.0
	236	ENDIF
	237
	238	ENDIF
	239
	240	!
[1019]	241	!-- Arrays needed for reasons of speed optimization for cache version.
	242	!-- For the vector version the buffer arrays are not necessary,
[713]	243	!-- because the the fluxes can swapped directly inside the loops of the
	244	!-- advection routines.
	245	IF ( loop_optimization /= 'vector' ) THEN
	246
	247	IF ( ws_scheme_mom ) THEN
	248
[736]	249	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
	250	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
	251	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
	252	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
	253	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
	254	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
	255	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	256	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	257	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	258	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	259	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	260	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
[713]	261
	262	ENDIF
	263
	264	IF ( ws_scheme_sca ) THEN
	265
[736]	266	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
	267	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
	268	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
	269	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
	270	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	271	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	272	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	273	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
[713]	274
[1115]	275	IF ( humidity .OR. passive_scalar ) THEN
[1053]	276	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
[736]	277	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
[1053]	278	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
[736]	279	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
[713]	280	ENDIF
	281
[1115]	282	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	283	precipitation ) THEN
[1053]	284	ALLOCATE( flux_s_qr(nzb+1:nzt,0:threads_per_task-1), &
	285	diss_s_qr(nzb+1:nzt,0:threads_per_task-1), &
	286	flux_s_nr(nzb+1:nzt,0:threads_per_task-1), &
	287	diss_s_nr(nzb+1:nzt,0:threads_per_task-1) )
	288	ALLOCATE( flux_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	289	diss_l_qr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	290	flux_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	291	diss_l_nr(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
[1115]	292	ENDIF
[1053]	293
[713]	294	IF ( ocean ) THEN
[736]	295	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), &
	296	diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
	297	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
	298	diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
[713]	299	ENDIF
	300
	301	ENDIF
	302
	303	ENDIF
	304
	305	END SUBROUTINE ws_init
	306
	307
	308	!------------------------------------------------------------------------------!
[1115]	309	! Initialize variables used for storing statistic quantities (fluxes, variances)
[713]	310	!------------------------------------------------------------------------------!
	311	SUBROUTINE ws_statistics
	312
	313	USE control_parameters
	314	USE statistics
	315
	316	IMPLICIT NONE
	317
	318	!
	319	!-- The arrays needed for statistical evaluation are set to to 0 at the
[861]	320	!-- beginning of prognostic_equations.
[713]	321	IF ( ws_scheme_mom ) THEN
	322	sums_wsus_ws_l = 0.0
	323	sums_wsvs_ws_l = 0.0
	324	sums_us2_ws_l = 0.0
	325	sums_vs2_ws_l = 0.0
	326	sums_ws2_ws_l = 0.0
	327	ENDIF
	328
	329	IF ( ws_scheme_sca ) THEN
	330	sums_wspts_ws_l = 0.0
[1115]	331	IF ( humidity .OR. passive_scalar ) sums_wsqs_ws_l = 0.0
	332	IF ( cloud_physics .AND. icloud_scheme == 0 .AND. &
	333	precipitation ) THEN
[1053]	334	sums_wsqrs_ws_l = 0.0
	335	sums_wsnrs_ws_l = 0.0
	336	ENDIF
[713]	337	IF ( ocean ) sums_wssas_ws_l = 0.0
	338
	339	ENDIF
	340
	341	END SUBROUTINE ws_statistics
	342
	343
	344	!------------------------------------------------------------------------------!
	345	! Scalar advection - Call for grid point i,j
	346	!------------------------------------------------------------------------------!
[1053]	347	SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char, swap_flux_y_local, &
	348	swap_diss_y_local, swap_flux_x_local, &
[736]	349	swap_diss_x_local, i_omp, tn )
[713]	350
	351	USE arrays_3d
	352	USE constants
	353	USE control_parameters
	354	USE grid_variables
	355	USE indices
[736]	356	USE pegrid
[713]	357	USE statistics
	358
	359	IMPLICIT NONE
	360
[888]	361	INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, &
	362	ibit7, ibit8, i_omp, j, k, k_mm, k_pp, k_ppp, tn
[862]	363	REAL :: diss_d, div, flux_d, u_comp, v_comp
[1010]	364	#if defined( __nopointer )
	365	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
	366	#else
	367	REAL, DIMENSION(:,:,:), POINTER :: sk
	368	#endif
	369	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, &
	370	flux_r, flux_t
	371	REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_diss_y_local, &
	372	swap_flux_y_local
	373	REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: &
	374	swap_diss_x_local, &
[862]	375	swap_flux_x_local
[1010]	376	CHARACTER (LEN = *), INTENT(IN) :: sk_char
[713]	377
[861]	378	!
	379	!-- Compute southside fluxes of the respective PE bounds.
	380	IF ( j == nys ) THEN
	381	!
	382	!-- Up to the top of the highest topography.
	383	DO k = nzb+1, nzb_max
[713]	384
[888]	385	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
	386	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
	387	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
	388
[861]	389	v_comp = v(k,j,i) - v_gtrans
	390	swap_flux_y_local(k,tn) = v_comp * ( &
[888]	391	( 37.0 * ibit5 * adv_sca_5 &
	392	+ 7.0 * ibit4 * adv_sca_3 &
	393	+ ibit3 * adv_sca_1 &
	394	) * &
	395	( sk(k,j,i) + sk(k,j-1,i) ) &
	396	- ( 8.0 * ibit5 * adv_sca_5 &
	397	+ ibit4 * adv_sca_3 &
	398	) * &
	399	( sk(k,j+1,i) + sk(k,j-2,i) ) &
	400	+ ( ibit5 * adv_sca_5 &
	401	) * &
	402	( sk(k,j+2,i) + sk(k,j-3,i) ) &
[861]	403	)
[713]	404
[861]	405	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
[888]	406	( 10.0 * ibit5 * adv_sca_5 &
	407	+ 3.0 * ibit4 * adv_sca_3 &
	408	+ ibit3 * adv_sca_1 &
	409	) * &
	410	( sk(k,j,i) - sk(k,j-1,i) ) &
	411	- ( 5.0 * ibit5 * adv_sca_5 &
	412	+ ibit4 * adv_sca_3 &
	413	) * &
	414	( sk(k,j+1,i) - sk(k,j-2,i) ) &
	415	+ ( ibit5 * adv_sca_5 &
	416	) * &
	417	( sk(k,j+2,i) - sk(k,j-3,i) ) &
[861]	418	)
[713]	419
[861]	420	ENDDO
	421	!
	422	!-- Above to the top of the highest topography. No degradation necessary.
	423	DO k = nzb_max+1, nzt
[713]	424
[888]	425	v_comp = v(k,j,i) - v_gtrans
[861]	426	swap_flux_y_local(k,tn) = v_comp * ( &
	427	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
	428	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
	429	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
	430	) * adv_sca_5
	431	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
	432	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
	433	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
	434	+ sk(k,j+2,i) - sk(k,j-3,i) &
	435	) * adv_sca_5
[713]	436
[861]	437	ENDDO
[713]	438
[861]	439	ENDIF
	440	!
	441	!-- Compute leftside fluxes of the respective PE bounds.
	442	IF ( i == i_omp ) THEN
	443
	444	DO k = nzb+1, nzb_max
[713]	445
[888]	446	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
	447	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
	448	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
	449
	450	u_comp = u(k,j,i) - u_gtrans
[861]	451	swap_flux_x_local(k,j,tn) = u_comp * ( &
[888]	452	( 37.0 * ibit2 * adv_sca_5 &
	453	+ 7.0 * ibit1 * adv_sca_3 &
	454	+ ibit0 * adv_sca_1 &
	455	) * &
	456	( sk(k,j,i) + sk(k,j,i-1) ) &
	457	- ( 8.0 * ibit2 * adv_sca_5 &
	458	+ ibit1 * adv_sca_3 &
	459	) * &
	460	( sk(k,j,i+1) + sk(k,j,i-2) ) &
	461	+ ( ibit2 * adv_sca_5 &
	462	) * &
	463	( sk(k,j,i+2) + sk(k,j,i-3) ) &
[861]	464	)
[713]	465
[861]	466	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
[888]	467	( 10.0 * ibit2 * adv_sca_5 &
	468	+ 3.0 * ibit1 * adv_sca_3 &
	469	+ ibit0 * adv_sca_1 &
	470	) * &
	471	( sk(k,j,i) - sk(k,j,i-1) ) &
	472	- ( 5.0 * ibit2 * adv_sca_5 &
	473	+ ibit1 * adv_sca_3 &
	474	) * &
	475	( sk(k,j,i+1) - sk(k,j,i-2) ) &
	476	+ ( ibit2 * adv_sca_5 &
	477	) * &
	478	( sk(k,j,i+2) - sk(k,j,i-3) ) &
[861]	479	)
[713]	480
[861]	481	ENDDO
[713]	482
[861]	483	DO k = nzb_max+1, nzt
[713]	484
[861]	485	u_comp = u(k,j,i) - u_gtrans
	486	swap_flux_x_local(k,j,tn) = u_comp * ( &
	487	37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) &
	488	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
	489	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
	490	) * adv_sca_5
[713]	491
[861]	492	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
	493	10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) &
	494	- 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
	495	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
	496	) * adv_sca_5
[713]	497
[861]	498	ENDDO
	499
	500	ENDIF
[713]	501
[861]	502	flux_t(0) = 0.0
	503	diss_t(0) = 0.0
	504	flux_d = 0.0
	505	diss_d = 0.0
	506	!
	507	!-- Now compute the fluxes and tendency terms for the horizontal and
	508	!-- vertical parts up to the top of the highest topography.
	509	DO k = nzb+1, nzb_max
	510	!
	511	!-- Note: It is faster to conduct all multiplications explicitly, e.g.
	512	!-- * adv_sca_5 ... than to determine a factor and multiplicate the
[888]	513	!-- flux at the end.
	514
	515	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
	516	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
	517	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
	518
[861]	519	u_comp = u(k,j,i+1) - u_gtrans
	520	flux_r(k) = u_comp * ( &
[888]	521	( 37.0 * ibit2 * adv_sca_5 &
	522	+ 7.0 * ibit1 * adv_sca_3 &
	523	+ ibit0 * adv_sca_1 &
[861]	524	) * &
	525	( sk(k,j,i+1) + sk(k,j,i) ) &
[888]	526	- ( 8.0 * ibit2 * adv_sca_5 &
	527	+ ibit1 * adv_sca_3 &
[861]	528	) * &
	529	( sk(k,j,i+2) + sk(k,j,i-1) ) &
[888]	530	+ ( ibit2 * adv_sca_5 &
	531	) * &
	532	( sk(k,j,i+3) + sk(k,j,i-2) ) &
[861]	533	)
[713]	534
[861]	535	diss_r(k) = -ABS( u_comp ) * ( &
[888]	536	( 10.0 * ibit2 * adv_sca_5 &
	537	+ 3.0 * ibit1 * adv_sca_3 &
	538	+ ibit0 * adv_sca_1 &
[861]	539	) * &
	540	( sk(k,j,i+1) - sk(k,j,i) ) &
[888]	541	- ( 5.0 * ibit2 * adv_sca_5 &
	542	+ ibit1 * adv_sca_3 &
[861]	543	) * &
	544	( sk(k,j,i+2) - sk(k,j,i-1) ) &
[888]	545	+ ( ibit2 * adv_sca_5 &
[861]	546	) * &
	547	( sk(k,j,i+3) - sk(k,j,i-2) ) &
	548	)
[713]	549
[888]	550	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
	551	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
	552	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
	553
[861]	554	v_comp = v(k,j+1,i) - v_gtrans
	555	flux_n(k) = v_comp * ( &
[888]	556	( 37.0 * ibit5 * adv_sca_5 &
	557	+ 7.0 * ibit4 * adv_sca_3 &
	558	+ ibit3 * adv_sca_1 &
[861]	559	) * &
	560	( sk(k,j+1,i) + sk(k,j,i) ) &
[888]	561	- ( 8.0 * ibit5 * adv_sca_5 &
	562	+ ibit4 * adv_sca_3 &
[861]	563	) * &
	564	( sk(k,j+2,i) + sk(k,j-1,i) ) &
[888]	565	+ ( ibit5 * adv_sca_5 &
[861]	566	) * &
	567	( sk(k,j+3,i) + sk(k,j-2,i) ) &
	568	)
[713]	569
[861]	570	diss_n(k) = -ABS( v_comp ) * ( &
[888]	571	( 10.0 * ibit5 * adv_sca_5 &
	572	+ 3.0 * ibit4 * adv_sca_3 &
	573	+ ibit3 * adv_sca_1 &
[861]	574	) * &
	575	( sk(k,j+1,i) - sk(k,j,i) ) &
[888]	576	- ( 5.0 * ibit5 * adv_sca_5 &
	577	+ ibit4 * adv_sca_3 &
[861]	578	) * &
	579	( sk(k,j+2,i) - sk(k,j-1,i) ) &
[888]	580	+ ( ibit5 * adv_sca_5 &
[861]	581	) * &
	582	( sk(k,j+3,i) - sk(k,j-2,i) ) &
	583	)
	584	!
	585	!-- k index has to be modified near bottom and top, else array
	586	!-- subscripts will be exceeded.
[888]	587	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
	588	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
	589	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
[713]	590
[888]	591	k_ppp = k + 3 * ibit8
	592	k_pp = k + 2 * ( 1 - ibit6 )
	593	k_mm = k - 2 * ibit8
[713]	594
[888]	595
[861]	596	flux_t(k) = w(k,j,i) * ( &
[888]	597	( 37.0 * ibit8 * adv_sca_5 &
	598	+ 7.0 * ibit7 * adv_sca_3 &
	599	+ ibit6 * adv_sca_1 &
[861]	600	) * &
	601	( sk(k+1,j,i) + sk(k,j,i) ) &
[888]	602	- ( 8.0 * ibit8 * adv_sca_5 &
	603	+ ibit7 * adv_sca_3 &
[861]	604	) * &
	605	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
[888]	606	+ ( ibit8 * adv_sca_5 &
[861]	607	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
	608	)
[713]	609
[861]	610	diss_t(k) = -ABS( w(k,j,i) ) * ( &
[888]	611	( 10.0 * ibit8 * adv_sca_5 &
	612	+ 3.0 * ibit7 * adv_sca_3 &
	613	+ ibit6 * adv_sca_1 &
[861]	614	) * &
	615	( sk(k+1,j,i) - sk(k,j,i) ) &
[888]	616	- ( 5.0 * ibit8 * adv_sca_5 &
	617	+ ibit7 * adv_sca_3 &
[861]	618	) * &
	619	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
[888]	620	+ ( ibit8 * adv_sca_5 &
[861]	621	) * &
	622	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
	623	)
[713]	624	!
[861]	625	!-- Calculate the divergence of the velocity field. A respective
	626	!-- correction is needed to overcome numerical instabilities caused
	627	!-- by a not sufficient reduction of divergences near topography.
	628	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
	629	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
	630	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[713]	631
[861]	632	tend(k,j,i) = tend(k,j,i) - ( &
	633	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
	634	swap_diss_x_local(k,j,tn) ) * ddx &
	635	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
	636	swap_diss_y_local(k,tn) ) * ddy &
	637	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
	638	) * ddzw(k) &
	639	) + sk(k,j,i) * div
[713]	640
[736]	641	swap_flux_y_local(k,tn) = flux_n(k)
	642	swap_diss_y_local(k,tn) = diss_n(k)
	643	swap_flux_x_local(k,j,tn) = flux_r(k)
	644	swap_diss_x_local(k,j,tn) = diss_r(k)
[861]	645	flux_d = flux_t(k)
	646	diss_d = diss_t(k)
	647
[713]	648	ENDDO
	649	!
[861]	650	!-- Now compute the fluxes and tendency terms for the horizontal and
	651	!-- vertical parts above the top of the highest topography. No degradation
	652	!-- for the horizontal parts, but for the vertical it is stell needed.
	653	DO k = nzb_max+1, nzt
[713]	654
[861]	655	u_comp = u(k,j,i+1) - u_gtrans
	656	flux_r(k) = u_comp * ( &
	657	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
	658	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
	659	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
	660	diss_r(k) = -ABS( u_comp ) * ( &
	661	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
	662	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
	663	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
[713]	664
[861]	665	v_comp = v(k,j+1,i) - v_gtrans
	666	flux_n(k) = v_comp * ( &
	667	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
	668	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
	669	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
	670	diss_n(k) = -ABS( v_comp ) * ( &
	671	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
	672	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
	673	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
[713]	674	!
[861]	675	!-- k index has to be modified near bottom and top, else array
	676	!-- subscripts will be exceeded.
[888]	677	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
	678	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
	679	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
[713]	680
[888]	681	k_ppp = k + 3 * ibit8
	682	k_pp = k + 2 * ( 1 - ibit6 )
	683	k_mm = k - 2 * ibit8
	684
	685
[861]	686	flux_t(k) = w(k,j,i) * ( &
[888]	687	( 37.0 * ibit8 * adv_sca_5 &
	688	+ 7.0 * ibit7 * adv_sca_3 &
	689	+ ibit6 * adv_sca_1 &
[861]	690	) * &
	691	( sk(k+1,j,i) + sk(k,j,i) ) &
[888]	692	- ( 8.0 * ibit8 * adv_sca_5 &
	693	+ ibit7 * adv_sca_3 &
[861]	694	) * &
	695	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
[888]	696	+ ( ibit8 * adv_sca_5 &
[861]	697	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
	698	)
	699
	700	diss_t(k) = -ABS( w(k,j,i) ) * ( &
[888]	701	( 10.0 * ibit8 * adv_sca_5 &
	702	+ 3.0 * ibit7 * adv_sca_3 &
	703	+ ibit6 * adv_sca_1 &
[861]	704	) * &
	705	( sk(k+1,j,i) - sk(k,j,i) ) &
[888]	706	- ( 5.0 * ibit8 * adv_sca_5 &
	707	+ ibit7 * adv_sca_3 &
[861]	708	) * &
	709	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
[888]	710	+ ( ibit8 * adv_sca_5 &
[861]	711	) * &
	712	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
	713	)
[713]	714	!
[861]	715	!-- Calculate the divergence of the velocity field. A respective
	716	!-- correction is needed to overcome numerical instabilities introduced
	717	!-- by a not sufficient reduction of divergences near topography.
	718	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
	719	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
	720	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[713]	721
[861]	722	tend(k,j,i) = tend(k,j,i) - ( &
	723	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
	724	swap_diss_x_local(k,j,tn) ) * ddx &
	725	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
	726	swap_diss_y_local(k,tn) ) * ddy &
	727	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
	728	) * ddzw(k) &
	729	) + sk(k,j,i) * div
[713]	730
[861]	731	swap_flux_y_local(k,tn) = flux_n(k)
	732	swap_diss_y_local(k,tn) = diss_n(k)
	733	swap_flux_x_local(k,j,tn) = flux_r(k)
	734	swap_diss_x_local(k,j,tn) = diss_r(k)
	735	flux_d = flux_t(k)
	736	diss_d = diss_t(k)
	737
[713]	738	ENDDO
	739
	740	!
	741	!-- Evaluation of statistics
	742	SELECT CASE ( sk_char )
	743
	744	CASE ( 'pt' )
	745
[861]	746	DO k = nzb, nzt
	747	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
[713]	748	( flux_t(k) + diss_t(k) ) &
[743]	749	* weight_substep(intermediate_timestep_count)
[713]	750	ENDDO
	751
	752	CASE ( 'sa' )
	753
[861]	754	DO k = nzb, nzt
	755	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
[713]	756	( flux_t(k) + diss_t(k) ) &
[743]	757	* weight_substep(intermediate_timestep_count)
[713]	758	ENDDO
	759
	760	CASE ( 'q' )
	761
[861]	762	DO k = nzb, nzt
	763	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
[713]	764	( flux_t(k) + diss_t(k) ) &
[743]	765	* weight_substep(intermediate_timestep_count)
[713]	766	ENDDO
	767
[1053]	768	CASE ( 'qr' )
	769
	770	DO k = nzb, nzt
	771	sums_wsqrs_ws_l(k,tn) = sums_wsqrs_ws_l(k,tn) + &
	772	( flux_t(k) + diss_t(k) ) &
	773	* weight_substep(intermediate_timestep_count)
	774	ENDDO
	775
	776	CASE ( 'nr' )
	777
	778	DO k = nzb, nzt
	779	sums_wsnrs_ws_l(k,tn) = sums_wsnrs_ws_l(k,tn) + &
	780	( flux_t(k) + diss_t(k) ) &
	781	* weight_substep(intermediate_timestep_count)
	782	ENDDO
	783
[713]	784	END SELECT
	785
	786	END SUBROUTINE advec_s_ws_ij
	787
	788
	789
	790
	791	!------------------------------------------------------------------------------!
	792	! Advection of u-component - Call for grid point i,j
	793	!------------------------------------------------------------------------------!
[736]	794	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
[713]	795
	796	USE arrays_3d
	797	USE constants
	798	USE control_parameters
	799	USE grid_variables
	800	USE indices
	801	USE statistics
	802
	803	IMPLICIT NONE
	804
[888]	805	INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, &
	806	ibit16, ibit17, i_omp, j, k, k_mm, k_pp, k_ppp, tn
[862]	807	REAL :: diss_d, div, flux_d, gu, gv, u_comp_l, v_comp, w_comp
	808	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
	809	flux_t, u_comp
[713]	810
	811	gu = 2.0 * u_gtrans
	812	gv = 2.0 * v_gtrans
[861]	813	!
	814	!-- Compute southside fluxes for the respective boundary of PE
	815	IF ( j == nys ) THEN
	816
	817	DO k = nzb+1, nzb_max
[713]	818
[888]	819	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
	820	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
	821	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
	822
[713]	823	v_comp = v(k,j,i) + v(k,j,i-1) - gv
[861]	824	flux_s_u(k,tn) = v_comp * ( &
[888]	825	( 37.0 * ibit14 * adv_mom_5 &
	826	+ 7.0 * ibit13 * adv_mom_3 &
	827	+ ibit12 * adv_mom_1 &
	828	) * &
[861]	829	( u(k,j,i) + u(k,j-1,i) ) &
[888]	830	- ( 8.0 * ibit14 * adv_mom_5 &
	831	+ ibit13 * adv_mom_3 &
	832	) * &
[861]	833	( u(k,j+1,i) + u(k,j-2,i) ) &
[888]	834	+ ( ibit14 * adv_mom_5 &
	835	) * &
[861]	836	( u(k,j+2,i) + u(k,j-3,i) ) &
	837	)
	838
	839	diss_s_u(k,tn) = - ABS ( v_comp ) * ( &
[888]	840	( 10.0 * ibit14 * adv_mom_5 &
	841	+ 3.0 * ibit13 * adv_mom_3 &
	842	+ ibit12 * adv_mom_1 &
	843	) * &
[861]	844	( u(k,j,i) - u(k,j-1,i) ) &
[888]	845	- ( 5.0 * ibit14 * adv_mom_5 &
	846	+ ibit13 * adv_mom_3 &
	847	) * &
[861]	848	( u(k,j+1,i) - u(k,j-2,i) ) &
[888]	849	+ ( ibit14 * adv_mom_5 &
	850	) * &
[861]	851	( u(k,j+2,i) - u(k,j-3,i) ) &
	852	)
	853
	854	ENDDO
	855
	856	DO k = nzb_max+1, nzt
	857
	858	v_comp = v(k,j,i) + v(k,j,i-1) - gv
	859	flux_s_u(k,tn) = v_comp * ( &
[713]	860	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
	861	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
	862	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
[861]	863	diss_s_u(k,tn) = - ABS(v_comp) * ( &
[713]	864	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
	865	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
	866	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
[861]	867
[713]	868	ENDDO
	869
	870	ENDIF
[861]	871	!
	872	!-- Compute leftside fluxes for the respective boundary of PE
	873	IF ( i == i_omp ) THEN
[713]	874
[861]	875	DO k = nzb+1, nzb_max
	876
[888]	877	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
	878	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
	879	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
	880
[861]	881	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
	882	flux_l_u(k,j,tn) = u_comp_l * ( &
[888]	883	( 37.0 * ibit11 * adv_mom_5 &
	884	+ 7.0 * ibit10 * adv_mom_3 &
	885	+ ibit9 * adv_mom_1 &
	886	) * &
[861]	887	( u(k,j,i) + u(k,j,i-1) ) &
[888]	888	- ( 8.0 * ibit11 * adv_mom_5 &
	889	+ ibit10 * adv_mom_3 &
	890	) * &
[861]	891	( u(k,j,i+1) + u(k,j,i-2) ) &
[888]	892	+ ( ibit11 * adv_mom_5 &
	893	) * &
[861]	894	( u(k,j,i+2) + u(k,j,i-3) ) &
	895	)
	896
	897	diss_l_u(k,j,tn) = - ABS( u_comp_l ) * ( &
[888]	898	( 10.0 * ibit11 * adv_mom_5 &
	899	+ 3.0 * ibit10 * adv_mom_3 &
	900	+ ibit9 * adv_mom_1 &
	901	) * &
[861]	902	( u(k,j,i) - u(k,j,i-1) ) &
[888]	903	- ( 5.0 * ibit11 * adv_mom_5 &
	904	+ ibit10 * adv_mom_3 &
	905	) * &
[861]	906	( u(k,j,i+1) - u(k,j,i-2) ) &
[888]	907	+ ( ibit11 * adv_mom_5 &
	908	) * &
[861]	909	( u(k,j,i+2) - u(k,j,i-3) ) &
	910	)
	911
	912	ENDDO
	913
	914	DO k = nzb_max+1, nzt
	915
	916	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
	917	flux_l_u(k,j,tn) = u_comp_l * ( &
	918	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
	919	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
[713]	920	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
[861]	921	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
	922	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
	923	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
[713]	924	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
[861]	925
[713]	926	ENDDO
	927
	928	ENDIF
[861]	929
	930	flux_t(0) = 0.0
	931	diss_t(0) = 0.0
	932	flux_d = 0.0
	933	diss_d = 0.0
	934	!
	935	!-- Now compute the fluxes tendency terms for the horizontal and
	936	!-- vertical parts.
	937	DO k = nzb+1, nzb_max
	938
[888]	939	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
	940	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
	941	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
	942
[861]	943	u_comp(k) = u(k,j,i+1) + u(k,j,i)
	944	flux_r(k) = ( u_comp(k) - gu ) * ( &
[888]	945	( 37.0 * ibit11 * adv_mom_5 &
	946	+ 7.0 * ibit10 * adv_mom_3 &
	947	+ ibit9 * adv_mom_1 &
[861]	948	) * &
	949	( u(k,j,i+1) + u(k,j,i) ) &
[888]	950	- ( 8.0 * ibit11 * adv_mom_5 &
	951	+ ibit10 * adv_mom_3 &
[861]	952	) * &
	953	( u(k,j,i+2) + u(k,j,i-1) ) &
[888]	954	+ ( ibit11 * adv_mom_5 &
[861]	955	) * &
	956	( u(k,j,i+3) + u(k,j,i-2) ) &
	957	)
	958
	959	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
[888]	960	( 10.0 * ibit11 * adv_mom_5 &
	961	+ 3.0 * ibit10 * adv_mom_3 &
	962	+ ibit9 * adv_mom_1 &
[861]	963	) * &
	964	( u(k,j,i+1) - u(k,j,i) ) &
[888]	965	- ( 5.0 * ibit11 * adv_mom_5 &
	966	+ ibit10 * adv_mom_3 &
[861]	967	) * &
	968	( u(k,j,i+2) - u(k,j,i-1) ) &
[888]	969	+ ( ibit11 * adv_mom_5 &
[861]	970	) * &
	971	( u(k,j,i+3) - u(k,j,i-2) ) &
	972	)
	973
[888]	974	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
	975	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
	976	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
	977
[861]	978	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
	979	flux_n(k) = v_comp * ( &
[888]	980	( 37.0 * ibit14 * adv_mom_5 &
	981	+ 7.0 * ibit13 * adv_mom_3 &
	982	+ ibit12 * adv_mom_1 &
[861]	983	) * &
	984	( u(k,j+1,i) + u(k,j,i) ) &
[888]	985	- ( 8.0 * ibit14 * adv_mom_5 &
	986	+ ibit13 * adv_mom_3 &
[861]	987	) * &
	988	( u(k,j+2,i) + u(k,j-1,i) ) &
[888]	989	+ ( ibit14 * adv_mom_5 &
[861]	990	) * &
	991	( u(k,j+3,i) + u(k,j-2,i) ) &
[888]	992	)
[861]	993
	994	diss_n(k) = - ABS ( v_comp ) * ( &
[888]	995	( 10.0 * ibit14 * adv_mom_5 &
	996	+ 3.0 * ibit13 * adv_mom_3 &
	997	+ ibit12 * adv_mom_1 &
[861]	998	) * &
	999	( u(k,j+1,i) - u(k,j,i) ) &
[888]	1000	- ( 5.0 * ibit14 * adv_mom_5 &
	1001	+ ibit13 * adv_mom_3 &
[861]	1002	) * &
	1003	( u(k,j+2,i) - u(k,j-1,i) ) &
[888]	1004	+ ( ibit14 * adv_mom_5 &
[861]	1005	) * &
	1006	( u(k,j+3,i) - u(k,j-2,i) ) &
[888]	1007	)
[861]	1008	!
	1009	!-- k index has to be modified near bottom and top, else array
	1010	!-- subscripts will be exceeded.
[888]	1011	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
	1012	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
	1013	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
[861]	1014
[888]	1015	k_ppp = k + 3 * ibit17
	1016	k_pp = k + 2 * ( 1 - ibit15 )
	1017	k_mm = k - 2 * ibit17
	1018
[861]	1019	w_comp = w(k,j,i) + w(k,j,i-1)
	1020	flux_t(k) = w_comp * ( &
[888]	1021	( 37.0 * ibit17 * adv_mom_5 &
	1022	+ 7.0 * ibit16 * adv_mom_3 &
	1023	+ ibit15 * adv_mom_1 &
[861]	1024	) * &
	1025	( u(k+1,j,i) + u(k,j,i) ) &
[888]	1026	- ( 8.0 * ibit17 * adv_mom_5 &
	1027	+ ibit16 * adv_mom_3 &
[861]	1028	) * &
	1029	( u(k_pp,j,i) + u(k-1,j,i) ) &
[888]	1030	+ ( ibit17 * adv_mom_5 &
[861]	1031	) * &
	1032	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
	1033	)
	1034
	1035	diss_t(k) = - ABS( w_comp ) * ( &
[888]	1036	( 10.0 * ibit17 * adv_mom_5 &
	1037	+ 3.0 * ibit16 * adv_mom_3 &
	1038	+ ibit15 * adv_mom_1 &
[861]	1039	) * &
	1040	( u(k+1,j,i) - u(k,j,i) ) &
[888]	1041	- ( 5.0 * ibit17 * adv_mom_5 &
	1042	+ ibit16 * adv_mom_3 &
[861]	1043	) * &
	1044	( u(k_pp,j,i) - u(k-1,j,i) ) &
[888]	1045	+ ( ibit17 * adv_mom_5 &
[861]	1046	) * &
	1047	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
	1048	)
	1049	!
	1050	!-- Calculate the divergence of the velocity field. A respective
	1051	!-- correction is needed to overcome numerical instabilities introduced
	1052	!-- by a not sufficient reduction of divergences near topography.
	1053	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
	1054	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
	1055	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) &
	1056	) * 0.5
	1057
[713]	1058	tend(k,j,i) = tend(k,j,i) - ( &
	1059	( flux_r(k) + diss_r(k) &
[861]	1060	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
[713]	1061	+ ( flux_n(k) + diss_n(k) &
[861]	1062	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
	1063	+ ( flux_t(k) + diss_t(k) &
	1064	- flux_d - diss_d ) * ddzw(k) &
	1065	) + div * u(k,j,i)
[713]	1066
[736]	1067	flux_l_u(k,j,tn) = flux_r(k)
	1068	diss_l_u(k,j,tn) = diss_r(k)
	1069	flux_s_u(k,tn) = flux_n(k)
	1070	diss_s_u(k,tn) = diss_n(k)
[861]	1071	flux_d = flux_t(k)
	1072	diss_d = diss_t(k)
[713]	1073	!
[861]	1074	!-- Statistical Evaluation of u'u'. The factor has to be applied for
[713]	1075	!-- right evaluation when gallilei_trans = .T. .
[861]	1076	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
	1077	+ ( flux_r(k) * &
	1078	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	1079	/ ( u_comp(k) - gu + 1.0E-20 ) &
	1080	+ diss_r(k) * &
	1081	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	1082	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
	1083	* weight_substep(intermediate_timestep_count)
	1084	!
	1085	!-- Statistical Evaluation of w'u'.
	1086	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
	1087	+ ( flux_t(k) + diss_t(k) ) &
	1088	* weight_substep(intermediate_timestep_count)
	1089	ENDDO
[713]	1090
[861]	1091	DO k = nzb_max+1, nzt
	1092
	1093	u_comp(k) = u(k,j,i+1) + u(k,j,i)
	1094	flux_r(k) = ( u_comp(k) - gu ) * ( &
	1095	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
	1096	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
	1097	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
	1098	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
	1099	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
	1100	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
	1101	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
	1102
	1103	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
	1104	flux_n(k) = v_comp * ( &
	1105	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
	1106	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
	1107	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
	1108	diss_n(k) = - ABS( v_comp ) * ( &
	1109	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
	1110	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
	1111	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
[713]	1112	!
[861]	1113	!-- k index has to be modified near bottom and top, else array
	1114	!-- subscripts will be exceeded.
[888]	1115	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
	1116	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
	1117	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
[713]	1118
[888]	1119	k_ppp = k + 3 * ibit17
	1120	k_pp = k + 2 * ( 1 - ibit15 )
	1121	k_mm = k - 2 * ibit17
	1122
[713]	1123	w_comp = w(k,j,i) + w(k,j,i-1)
[861]	1124	flux_t(k) = w_comp * ( &
[888]	1125	( 37.0 * ibit17 * adv_mom_5 &
	1126	+ 7.0 * ibit16 * adv_mom_3 &
	1127	+ ibit15 * adv_mom_1 &
[861]	1128	) * &
	1129	( u(k+1,j,i) + u(k,j,i) ) &
[888]	1130	- ( 8.0 * ibit17 * adv_mom_5 &
	1131	+ ibit16 * adv_mom_3 &
[861]	1132	) * &
	1133	( u(k_pp,j,i) + u(k-1,j,i) ) &
[888]	1134	+ ( ibit17 * adv_mom_5 &
[861]	1135	) * &
	1136	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
	1137	)
[713]	1138
[861]	1139	diss_t(k) = - ABS( w_comp ) * ( &
[888]	1140	( 10.0 * ibit17 * adv_mom_5 &
	1141	+ 3.0 * ibit16 * adv_mom_3 &
	1142	+ ibit15 * adv_mom_1 &
[861]	1143	) * &
	1144	( u(k+1,j,i) - u(k,j,i) ) &
[888]	1145	- ( 5.0 * ibit17 * adv_mom_5 &
	1146	+ ibit16 * adv_mom_3 &
[861]	1147	) * &
	1148	( u(k_pp,j,i) - u(k-1,j,i) ) &
[888]	1149	+ ( ibit17 * adv_mom_5 &
[861]	1150	) * &
	1151	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
	1152	)
[713]	1153	!
[861]	1154	!-- Calculate the divergence of the velocity field. A respective
	1155	!-- correction is needed to overcome numerical instabilities introduced
	1156	!-- by a not sufficient reduction of divergences near topography.
	1157	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
	1158	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
	1159	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) * ddzw(k) &
	1160	) * 0.5
[713]	1161
[861]	1162	tend(k,j,i) = tend(k,j,i) - ( &
	1163	( flux_r(k) + diss_r(k) &
	1164	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
	1165	+ ( flux_n(k) + diss_n(k) &
	1166	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &
	1167	+ ( flux_t(k) + diss_t(k) &
	1168	- flux_d - diss_d ) * ddzw(k) &
	1169	) + div * u(k,j,i)
[713]	1170
[861]	1171	flux_l_u(k,j,tn) = flux_r(k)
	1172	diss_l_u(k,j,tn) = diss_r(k)
	1173	flux_s_u(k,tn) = flux_n(k)
	1174	diss_s_u(k,tn) = diss_n(k)
	1175	flux_d = flux_t(k)
	1176	diss_d = diss_t(k)
[713]	1177	!
[861]	1178	!-- Statistical Evaluation of u'u'. The factor has to be applied for
	1179	!-- right evaluation when gallilei_trans = .T. .
	1180	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
	1181	+ ( flux_r(k) * &
	1182	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	1183	/ ( u_comp(k) - gu + 1.0E-20 ) &
	1184	+ diss_r(k) * &
	1185	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	1186	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
	1187	* weight_substep(intermediate_timestep_count)
	1188	!
	1189	!-- Statistical Evaluation of w'u'.
	1190	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
	1191	+ ( flux_t(k) + diss_t(k) ) &
	1192	* weight_substep(intermediate_timestep_count)
[713]	1193	ENDDO
	1194
[861]	1195	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
[713]	1196
[861]	1197
	1198
[713]	1199	END SUBROUTINE advec_u_ws_ij
	1200
	1201
	1202
[736]	1203	!-----------------------------------------------------------------------------!
[713]	1204	! Advection of v-component - Call for grid point i,j
[736]	1205	!-----------------------------------------------------------------------------!
	1206	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
[713]	1207
	1208	USE arrays_3d
	1209	USE constants
	1210	USE control_parameters
	1211	USE grid_variables
	1212	USE indices
	1213	USE statistics
	1214
	1215	IMPLICIT NONE
	1216
[888]	1217	INTEGER :: i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
	1218	ibit25, ibit26, i_omp, j, k, k_mm, k_pp, k_ppp, tn
[862]	1219	REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp_l, w_comp
	1220	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
	1221	flux_t, v_comp
[861]	1222
[713]	1223	gu = 2.0 * u_gtrans
	1224	gv = 2.0 * v_gtrans
	1225
	1226	!
[861]	1227	!-- Compute leftside fluxes for the respective boundary.
	1228	IF ( i == i_omp ) THEN
[736]	1229
[861]	1230	DO k = nzb+1, nzb_max
	1231
[888]	1232	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
	1233	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
	1234	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
	1235
[861]	1236	u_comp = u(k,j-1,i) + u(k,j,i) - gu
	1237	flux_l_v(k,j,tn) = u_comp * ( &
[888]	1238	( 37.0 * ibit20 * adv_mom_5 &
	1239	+ 7.0 * ibit19 * adv_mom_3 &
	1240	+ ibit18 * adv_mom_1 &
	1241	) * &
[861]	1242	( v(k,j,i) + v(k,j,i-1) ) &
[888]	1243	- ( 8.0 * ibit20 * adv_mom_5 &
	1244	+ ibit19 * adv_mom_3 &
	1245	) * &
[861]	1246	( v(k,j,i+1) + v(k,j,i-2) ) &
[888]	1247	+ ( ibit20 * adv_mom_5 &
	1248	) * &
[861]	1249	( v(k,j,i+2) + v(k,j,i-3) ) &
	1250	)
	1251
	1252	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
[888]	1253	( 10.0 * ibit20 * adv_mom_5 &
	1254	+ 3.0 * ibit19 * adv_mom_3 &
	1255	+ ibit18 * adv_mom_1 &
	1256	) * &
[861]	1257	( v(k,j,i) - v(k,j,i-1) ) &
[888]	1258	- ( 5.0 * ibit20 * adv_mom_5 &
	1259	+ ibit19 * adv_mom_3 &
	1260	) * &
[861]	1261	( v(k,j,i+1) - v(k,j,i-2) ) &
[888]	1262	+ ( ibit20 * adv_mom_5 &
	1263	) * &
[861]	1264	( v(k,j,i+2) - v(k,j,i-3) ) &
	1265	)
	1266
	1267	ENDDO
	1268
	1269	DO k = nzb_max+1, nzt
	1270
	1271	u_comp = u(k,j-1,i) + u(k,j,i) - gu
	1272	flux_l_v(k,j,tn) = u_comp * ( &
[713]	1273	37.0 * ( v(k,j,i) + v(k,j,i-1) ) &
	1274	- 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) &
	1275	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
[861]	1276	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
[713]	1277	10.0 * ( v(k,j,i) - v(k,j,i-1) ) &
	1278	- 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) &
	1279	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
[861]	1280
[713]	1281	ENDDO
	1282
	1283	ENDIF
[861]	1284	!
	1285	!-- Compute southside fluxes for the respective boundary.
	1286	IF ( j == nysv ) THEN
[713]	1287
[861]	1288	DO k = nzb+1, nzb_max
	1289
[888]	1290	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
	1291	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
	1292	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
	1293
[861]	1294	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
	1295	flux_s_v(k,tn) = v_comp_l * ( &
[888]	1296	( 37.0 * ibit23 * adv_mom_5 &
	1297	+ 7.0 * ibit22 * adv_mom_3 &
	1298	+ ibit21 * adv_mom_1 &
	1299	) * &
[861]	1300	( v(k,j,i) + v(k,j-1,i) ) &
[888]	1301	- ( 8.0 * ibit23 * adv_mom_5 &
	1302	+ ibit22 * adv_mom_3 &
	1303	) * &
[861]	1304	( v(k,j+1,i) + v(k,j-2,i) ) &
[888]	1305	+ ( ibit23 * adv_mom_5 &
	1306	) * &
[861]	1307	( v(k,j+2,i) + v(k,j-3,i) ) &
	1308	)
	1309
	1310	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
[888]	1311	( 10.0 * ibit23 * adv_mom_5 &
	1312	+ 3.0 * ibit22 * adv_mom_3 &
	1313	+ ibit21 * adv_mom_1 &
	1314	) * &
[861]	1315	( v(k,j,i) - v(k,j-1,i) ) &
[888]	1316	- ( 5.0 * ibit23 * adv_mom_5 &
	1317	+ ibit22 * adv_mom_3 &
	1318	) * &
[861]	1319	( v(k,j+1,i) - v(k,j-2,i) ) &
[888]	1320	+ ( ibit23 * adv_mom_5 &
	1321	) * &
[861]	1322	( v(k,j+2,i) - v(k,j-3,i) ) &
	1323	)
	1324
	1325	ENDDO
	1326
	1327	DO k = nzb_max+1, nzt
	1328
	1329	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
	1330	flux_s_v(k,tn) = v_comp_l * ( &
[713]	1331	37.0 * ( v(k,j,i) + v(k,j-1,i) ) &
	1332	- 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) &
	1333	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
[861]	1334	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
[713]	1335	10.0 * ( v(k,j,i) - v(k,j-1,i) ) &
	1336	- 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) &
	1337	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
[861]	1338
[713]	1339	ENDDO
	1340
	1341	ENDIF
[861]	1342
	1343	flux_t(0) = 0.0
	1344	diss_t(0) = 0.0
	1345	flux_d = 0.0
	1346	diss_d = 0.0
	1347	!
	1348	!-- Now compute the fluxes and tendency terms for the horizontal and
	1349	!-- verical parts.
	1350	DO k = nzb+1, nzb_max
	1351
[888]	1352	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
	1353	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
	1354	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
	1355
[861]	1356	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
	1357	flux_r(k) = u_comp * ( &
[888]	1358	( 37.0 * ibit20 * adv_mom_5 &
	1359	+ 7.0 * ibit19 * adv_mom_3 &
	1360	+ ibit18 * adv_mom_1 &
[861]	1361	) * &
	1362	( v(k,j,i+1) + v(k,j,i) ) &
[888]	1363	- ( 8.0 * ibit20 * adv_mom_5 &
	1364	+ ibit19 * adv_mom_3 &
[861]	1365	) * &
	1366	( v(k,j,i+2) + v(k,j,i-1) ) &
[888]	1367	+ ( ibit20 * adv_mom_5 &
[861]	1368	) * &
	1369	( v(k,j,i+3) + v(k,j,i-2) ) &
	1370	)
	1371
	1372	diss_r(k) = - ABS( u_comp ) * ( &
[888]	1373	( 10.0 * ibit20 * adv_mom_5 &
	1374	+ 3.0 * ibit19 * adv_mom_3 &
	1375	+ ibit18 * adv_mom_1 &
[861]	1376	) * &
	1377	( v(k,j,i+1) - v(k,j,i) ) &
[888]	1378	- ( 5.0 * ibit20 * adv_mom_5 &
	1379	+ ibit19 * adv_mom_3 &
[861]	1380	) * &
	1381	( v(k,j,i+2) - v(k,j,i-1) ) &
[888]	1382	+ ( ibit20 * adv_mom_5 &
[861]	1383	) * &
	1384	( v(k,j,i+3) - v(k,j,i-2) ) &
	1385	)
	1386
[888]	1387	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
	1388	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
	1389	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
	1390
	1391
[861]	1392	v_comp(k) = v(k,j+1,i) + v(k,j,i)
	1393	flux_n(k) = ( v_comp(k) - gv ) * ( &
[888]	1394	( 37.0 * ibit23 * adv_mom_5 &
	1395	+ 7.0 * ibit22 * adv_mom_3 &
	1396	+ ibit21 * adv_mom_1 &
[861]	1397	) * &
	1398	( v(k,j+1,i) + v(k,j,i) ) &
[888]	1399	- ( 8.0 * ibit23 * adv_mom_5 &
	1400	+ ibit22 * adv_mom_3 &
[861]	1401	) * &
	1402	( v(k,j+2,i) + v(k,j-1,i) ) &
[888]	1403	+ ( ibit23 * adv_mom_5 &
[861]	1404	) * &
	1405	( v(k,j+3,i) + v(k,j-2,i) ) &
	1406	)
	1407
	1408	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
[888]	1409	( 10.0 * ibit23 * adv_mom_5 &
	1410	+ 3.0 * ibit22 * adv_mom_3 &
	1411	+ ibit21 * adv_mom_1 &
[861]	1412	) * &
	1413	( v(k,j+1,i) - v(k,j,i) ) &
[888]	1414	- ( 5.0 * ibit23 * adv_mom_5 &
	1415	+ ibit22 * adv_mom_3 &
[861]	1416	) * &
	1417	( v(k,j+2,i) - v(k,j-1,i) ) &
[888]	1418	+ ( ibit23 * adv_mom_5 &
[861]	1419	) * &
	1420	( v(k,j+3,i) - v(k,j-2,i) ) &
	1421	)
	1422	!
	1423	!-- k index has to be modified near bottom and top, else array
	1424	!-- subscripts will be exceeded.
[888]	1425	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
	1426	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
	1427	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
[861]	1428
[888]	1429	k_ppp = k + 3 * ibit26
	1430	k_pp = k + 2 * ( 1 - ibit24 )
	1431	k_mm = k - 2 * ibit26
	1432
[861]	1433	w_comp = w(k,j-1,i) + w(k,j,i)
	1434	flux_t(k) = w_comp * ( &
[888]	1435	( 37.0 * ibit26 * adv_mom_5 &
	1436	+ 7.0 * ibit25 * adv_mom_3 &
	1437	+ ibit24 * adv_mom_1 &
[861]	1438	) * &
	1439	( v(k+1,j,i) + v(k,j,i) ) &
[888]	1440	- ( 8.0 * ibit26 * adv_mom_5 &
	1441	+ ibit25 * adv_mom_3 &
[861]	1442	) * &
	1443	( v(k_pp,j,i) + v(k-1,j,i) ) &
[888]	1444	+ ( ibit26 * adv_mom_5 &
[861]	1445	) * &
	1446	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
	1447	)
	1448
	1449	diss_t(k) = - ABS( w_comp ) * ( &
[888]	1450	( 10.0 * ibit26 * adv_mom_5 &
	1451	+ 3.0 * ibit25 * adv_mom_3 &
	1452	+ ibit24 * adv_mom_1 &
[861]	1453	) * &
	1454	( v(k+1,j,i) - v(k,j,i) ) &
[888]	1455	- ( 5.0 * ibit26 * adv_mom_5 &
	1456	+ ibit25 * adv_mom_3 &
[861]	1457	) * &
	1458	( v(k_pp,j,i) - v(k-1,j,i) ) &
[888]	1459	+ ( ibit26 * adv_mom_5 &
[861]	1460	) * &
	1461	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
	1462	)
	1463	!
	1464	!-- Calculate the divergence of the velocity field. A respective
	1465	!-- correction is needed to overcome numerical instabilities introduced
	1466	!-- by a not sufficient reduction of divergences near topography.
	1467	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
	1468	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
	1469	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) &
	1470	) * 0.5
	1471
[713]	1472	tend(k,j,i) = tend(k,j,i) - ( &
	1473	( flux_r(k) + diss_r(k) &
[861]	1474	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
[713]	1475	+ ( flux_n(k) + diss_n(k) &
[861]	1476	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
	1477	+ ( flux_t(k) + diss_t(k) &
	1478	- flux_d - diss_d ) * ddzw(k) &
	1479	) + v(k,j,i) * div
	1480
[736]	1481	flux_l_v(k,j,tn) = flux_r(k)
	1482	diss_l_v(k,j,tn) = diss_r(k)
	1483	flux_s_v(k,tn) = flux_n(k)
	1484	diss_s_v(k,tn) = diss_n(k)
[861]	1485	flux_d = flux_t(k)
	1486	diss_d = diss_t(k)
[713]	1487
	1488	!
[861]	1489	!-- Statistical Evaluation of v'v'. The factor has to be applied for
[713]	1490	!-- right evaluation when gallilei_trans = .T. .
[861]	1491	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
[713]	1492	+ ( flux_n(k) &
[743]	1493	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
[713]	1494	/ ( v_comp(k) - gv + 1.0E-20 ) &
	1495	+ diss_n(k) &
[743]	1496	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
[713]	1497	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
[743]	1498	* weight_substep(intermediate_timestep_count)
[861]	1499	!
	1500	!-- Statistical Evaluation of w'v'.
	1501	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
	1502	+ ( flux_t(k) + diss_t(k) ) &
	1503	* weight_substep(intermediate_timestep_count)
[713]	1504
	1505	ENDDO
	1506
[861]	1507	DO k = nzb_max+1, nzt
	1508
	1509	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
	1510	flux_r(k) = u_comp * ( &
	1511	37.0 * ( v(k,j,i+1) + v(k,j,i) ) &
	1512	- 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) &
	1513	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
	1514
	1515	diss_r(k) = - ABS( u_comp ) * ( &
	1516	10.0 * ( v(k,j,i+1) - v(k,j,i) ) &
	1517	- 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) &
	1518	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
	1519
	1520
	1521	v_comp(k) = v(k,j+1,i) + v(k,j,i)
	1522	flux_n(k) = ( v_comp(k) - gv ) * ( &
	1523	37.0 * ( v(k,j+1,i) + v(k,j,i) ) &
	1524	- 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) &
	1525	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
	1526
	1527	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
	1528	10.0 * ( v(k,j+1,i) - v(k,j,i) ) &
	1529	- 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) &
	1530	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
[713]	1531	!
[861]	1532	!-- k index has to be modified near bottom and top, else array
	1533	!-- subscripts will be exceeded.
[888]	1534	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
	1535	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
	1536	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
[713]	1537
[888]	1538	k_ppp = k + 3 * ibit26
	1539	k_pp = k + 2 * ( 1 - ibit24 )
	1540	k_mm = k - 2 * ibit26
	1541
[713]	1542	w_comp = w(k,j-1,i) + w(k,j,i)
[861]	1543	flux_t(k) = w_comp * ( &
[888]	1544	( 37.0 * ibit26 * adv_mom_5 &
	1545	+ 7.0 * ibit25 * adv_mom_3 &
	1546	+ ibit24 * adv_mom_1 &
[861]	1547	) * &
	1548	( v(k+1,j,i) + v(k,j,i) ) &
[888]	1549	- ( 8.0 * ibit26 * adv_mom_5 &
	1550	+ ibit25 * adv_mom_3 &
[861]	1551	) * &
	1552	( v(k_pp,j,i) + v(k-1,j,i) ) &
[888]	1553	+ ( ibit26 * adv_mom_5 &
[861]	1554	) * &
	1555	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
	1556	)
[713]	1557
[861]	1558	diss_t(k) = - ABS( w_comp ) * ( &
[888]	1559	( 10.0 * ibit26 * adv_mom_5 &
	1560	+ 3.0 * ibit25 * adv_mom_3 &
	1561	+ ibit24 * adv_mom_1 &
[861]	1562	) * &
	1563	( v(k+1,j,i) - v(k,j,i) ) &
[888]	1564	- ( 5.0 * ibit26 * adv_mom_5 &
	1565	+ ibit25 * adv_mom_3 &
[861]	1566	) * &
	1567	( v(k_pp,j,i) - v(k-1,j,i) ) &
[888]	1568	+ ( ibit26 * adv_mom_5 &
[861]	1569	) * &
	1570	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
	1571	)
[713]	1572	!
[861]	1573	!-- Calculate the divergence of the velocity field. A respective
	1574	!-- correction is needed to overcome numerical instabilities introduced
	1575	!-- by a not sufficient reduction of divergences near topography.
	1576	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
	1577	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
	1578	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) * ddzw(k) &
	1579	) * 0.5
	1580
	1581	tend(k,j,i) = tend(k,j,i) - ( &
	1582	( flux_r(k) + diss_r(k) &
	1583	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
	1584	+ ( flux_n(k) + diss_n(k) &
	1585	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &
	1586	+ ( flux_t(k) + diss_t(k) &
	1587	- flux_d - diss_d ) * ddzw(k) &
	1588	) + v(k,j,i) * div
	1589
	1590	flux_l_v(k,j,tn) = flux_r(k)
	1591	diss_l_v(k,j,tn) = diss_r(k)
	1592	flux_s_v(k,tn) = flux_n(k)
	1593	diss_s_v(k,tn) = diss_n(k)
	1594	flux_d = flux_t(k)
	1595	diss_d = diss_t(k)
	1596
[713]	1597	!
[861]	1598	!-- Statistical Evaluation of v'v'. The factor has to be applied for
	1599	!-- right evaluation when gallilei_trans = .T. .
	1600	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
	1601	+ ( flux_n(k) &
	1602	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
	1603	/ ( v_comp(k) - gv + 1.0E-20 ) &
	1604	+ diss_n(k) &
	1605	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
	1606	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
	1607	* weight_substep(intermediate_timestep_count)
	1608	!
	1609	!-- Statistical Evaluation of w'v'.
	1610	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
	1611	+ ( flux_t(k) + diss_t(k) ) &
	1612	* weight_substep(intermediate_timestep_count)
	1613
[713]	1614	ENDDO
[861]	1615	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
[713]	1616
[861]	1617
[713]	1618	END SUBROUTINE advec_v_ws_ij
	1619
	1620
	1621
	1622	!------------------------------------------------------------------------------!
	1623	! Advection of w-component - Call for grid point i,j
	1624	!------------------------------------------------------------------------------!
[736]	1625	SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
[713]	1626
	1627	USE arrays_3d
	1628	USE constants
	1629	USE control_parameters
	1630	USE grid_variables
	1631	USE indices
	1632	USE statistics
	1633
	1634	IMPLICIT NONE
	1635
[888]	1636	INTEGER :: i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
	1637	ibit34, ibit35, i_omp, j, k, k_mm, k_pp, k_ppp, tn
[862]	1638	REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp
	1639	REAL, DIMENSION(nzb:nzt+1) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
	1640	flux_t
[713]	1641
	1642	gu = 2.0 * u_gtrans
	1643	gv = 2.0 * v_gtrans
[736]	1644
[861]	1645	!
	1646	!-- Compute southside fluxes for the respective boundary.
	1647	IF ( j == nys ) THEN
	1648
	1649	DO k = nzb+1, nzb_max
[1221]	1650	ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
[888]	1651	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
	1652	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
[861]	1653
[862]	1654	v_comp = v(k+1,j,i) + v(k,j,i) - gv
[861]	1655	flux_s_w(k,tn) = v_comp * ( &
[888]	1656	( 37.0 * ibit32 * adv_mom_5 &
	1657	+ 7.0 * ibit31 * adv_mom_3 &
	1658	+ ibit30 * adv_mom_1 &
	1659	) * &
[861]	1660	( w(k,j,i) + w(k,j-1,i) ) &
[888]	1661	- ( 8.0 * ibit32 * adv_mom_5 &
	1662	+ ibit31 * adv_mom_3 &
	1663	) * &
[861]	1664	( w(k,j+1,i) + w(k,j-2,i) ) &
[888]	1665	+ ( ibit32 * adv_mom_5 &
	1666	) * &
[861]	1667	( w(k,j+2,i) + w(k,j-3,i) ) &
	1668	)
	1669
	1670	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
[888]	1671	( 10.0 * ibit32 * adv_mom_5 &
	1672	+ 3.0 * ibit31 * adv_mom_3 &
	1673	+ ibit30 * adv_mom_1 &
	1674	) * &
[861]	1675	( w(k,j,i) - w(k,j-1,i) ) &
[888]	1676	- ( 5.0 * ibit32 * adv_mom_5 &
	1677	+ ibit31 * adv_mom_3 &
	1678	) * &
[861]	1679	( w(k,j+1,i) - w(k,j-2,i) ) &
[888]	1680	+ ( ibit32 * adv_mom_5 &
	1681	) * &
[861]	1682	( w(k,j+2,i) - w(k,j-3,i) ) &
	1683	)
	1684
	1685	ENDDO
	1686
	1687	DO k = nzb_max+1, nzt
	1688
	1689	v_comp = v(k+1,j,i) + v(k,j,i) - gv
	1690	flux_s_w(k,tn) = v_comp * ( &
[713]	1691	37.0 * ( w(k,j,i) + w(k,j-1,i) ) &
	1692	- 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) &
	1693	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
[861]	1694	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
[713]	1695	10.0 * ( w(k,j,i) - w(k,j-1,i) ) &
	1696	- 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) &
	1697	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
[861]	1698
[713]	1699	ENDDO
[861]	1700
[713]	1701	ENDIF
[861]	1702	!
	1703	!-- Compute leftside fluxes for the respective boundary.
	1704	IF ( i == i_omp ) THEN
	1705
	1706	DO k = nzb+1, nzb_max
	1707
[888]	1708	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
	1709	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
	1710	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
	1711
[862]	1712	u_comp = u(k+1,j,i) + u(k,j,i) - gu
[861]	1713	flux_l_w(k,j,tn) = u_comp * ( &
[888]	1714	( 37.0 * ibit29 * adv_mom_5 &
	1715	+ 7.0 * ibit28 * adv_mom_3 &
	1716	+ ibit27 * adv_mom_1 &
	1717	) * &
[861]	1718	( w(k,j,i) + w(k,j,i-1) ) &
[888]	1719	- ( 8.0 * ibit29 * adv_mom_5 &
	1720	+ ibit28 * adv_mom_3 &
	1721	) * &
[861]	1722	( w(k,j,i+1) + w(k,j,i-2) ) &
[888]	1723	+ ( ibit29 * adv_mom_5 &
	1724	) * &
[861]	1725	( w(k,j,i+2) + w(k,j,i-3) ) &
	1726	)
	1727
	1728	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
[888]	1729	( 10.0 * ibit29 * adv_mom_5 &
	1730	+ 3.0 * ibit28 * adv_mom_3 &
	1731	+ ibit27 * adv_mom_1 &
	1732	) * &
[861]	1733	( w(k,j,i) - w(k,j,i-1) ) &
[888]	1734	- ( 5.0 * ibit29 * adv_mom_5 &
	1735	+ ibit28 * adv_mom_3 &
	1736	) * &
[861]	1737	( w(k,j,i+1) - w(k,j,i-2) ) &
[888]	1738	+ ( ibit29 * adv_mom_5 &
	1739	) * &
[861]	1740	( w(k,j,i+2) - w(k,j,i-3) ) &
	1741	)
	1742
	1743	ENDDO
	1744
	1745	DO k = nzb_max+1, nzt
	1746
	1747	u_comp = u(k+1,j,i) + u(k,j,i) - gu
	1748	flux_l_w(k,j,tn) = u_comp * ( &
[713]	1749	37.0 * ( w(k,j,i) + w(k,j,i-1) ) &
	1750	- 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) &
	1751	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
[861]	1752	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
[713]	1753	10.0 * ( w(k,j,i) - w(k,j,i-1) ) &
	1754	- 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) &
	1755	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
[861]	1756
[713]	1757	ENDDO
[861]	1758
[713]	1759	ENDIF
[861]	1760	!
[862]	1761	!-- The lower flux has to be calculated explicetely for the tendency at
	1762	!-- the first w-level. For topography wall this is done implicitely by
	1763	!-- wall_flags_0.
	1764	k = nzb + 1
	1765	w_comp = w(k,j,i) + w(k-1,j,i)
	1766	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
	1767	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
	1768	flux_d = flux_t(0)
	1769	diss_d = diss_t(0)
	1770	!
[861]	1771	!-- Now compute the fluxes and tendency terms for the horizontal
	1772	!-- and vertical parts.
	1773	DO k = nzb+1, nzb_max
	1774
[888]	1775	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
	1776	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
	1777	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
	1778
[861]	1779	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
	1780	flux_r(k) = u_comp * ( &
[888]	1781	( 37.0 * ibit29 * adv_mom_5 &
	1782	+ 7.0 * ibit28 * adv_mom_3 &
	1783	+ ibit27 * adv_mom_1 &
[861]	1784	) * &
	1785	( w(k,j,i+1) + w(k,j,i) ) &
[888]	1786	- ( 8.0 * ibit29 * adv_mom_5 &
	1787	+ ibit28 * adv_mom_3 &
[861]	1788	) * &
	1789	( w(k,j,i+2) + w(k,j,i-1) ) &
[888]	1790	+ ( ibit29 * adv_mom_5 &
[861]	1791	) * &
	1792	( w(k,j,i+3) + w(k,j,i-2) ) &
[888]	1793	)
[861]	1794
	1795	diss_r(k) = - ABS( u_comp ) * ( &
[888]	1796	( 10.0 * ibit29 * adv_mom_5 &
	1797	+ 3.0 * ibit28 * adv_mom_3 &
	1798	+ ibit27 * adv_mom_1 &
[861]	1799	) * &
	1800	( w(k,j,i+1) - w(k,j,i) ) &
[888]	1801	- ( 5.0 * ibit29 * adv_mom_5 &
	1802	+ ibit28 * adv_mom_3 &
[861]	1803	) * &
	1804	( w(k,j,i+2) - w(k,j,i-1) ) &
[888]	1805	+ ( ibit29 * adv_mom_5 &
[861]	1806	) * &
	1807	( w(k,j,i+3) - w(k,j,i-2) ) &
	1808	)
	1809
[1221]	1810	ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
[888]	1811	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
	1812	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
	1813
[861]	1814	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
	1815	flux_n(k) = v_comp * ( &
[888]	1816	( 37.0 * ibit32 * adv_mom_5 &
	1817	+ 7.0 * ibit31 * adv_mom_3 &
	1818	+ ibit30 * adv_mom_1 &
[861]	1819	) * &
	1820	( w(k,j+1,i) + w(k,j,i) ) &
[888]	1821	- ( 8.0 * ibit32 * adv_mom_5 &
	1822	+ ibit31 * adv_mom_3 &
[861]	1823	) * &
	1824	( w(k,j+2,i) + w(k,j-1,i) ) &
[888]	1825	+ ( ibit32 * adv_mom_5 &
[861]	1826	) * &
	1827	( w(k,j+3,i) + w(k,j-2,i) ) &
	1828	)
	1829
	1830	diss_n(k) = - ABS( v_comp ) * ( &
[888]	1831	( 10.0 * ibit32 * adv_mom_5 &
	1832	+ 3.0 * ibit31 * adv_mom_3 &
	1833	+ ibit30 * adv_mom_1 &
[861]	1834	) * &
	1835	( w(k,j+1,i) - w(k,j,i) ) &
[888]	1836	- ( 5.0 * ibit32 * adv_mom_5 &
	1837	+ ibit31 * adv_mom_3 &
[861]	1838	) * &
	1839	( w(k,j+2,i) - w(k,j-1,i) ) &
[888]	1840	+ ( ibit32 * adv_mom_5 &
[861]	1841	) * &
	1842	( w(k,j+3,i) - w(k,j-2,i) ) &
	1843	)
	1844	!
	1845	!-- k index has to be modified near bottom and top, else array
	1846	!-- subscripts will be exceeded.
[1221]	1847	ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
	1848	ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
	1849	ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
[861]	1850
[888]	1851	k_ppp = k + 3 * ibit35
	1852	k_pp = k + 2 * ( 1 - ibit33 )
	1853	k_mm = k - 2 * ibit35
	1854
[861]	1855	w_comp = w(k+1,j,i) + w(k,j,i)
	1856	flux_t(k) = w_comp * ( &
[888]	1857	( 37.0 * ibit35 * adv_mom_5 &
	1858	+ 7.0 * ibit34 * adv_mom_3 &
	1859	+ ibit33 * adv_mom_1 &
[861]	1860	) * &
	1861	( w(k+1,j,i) + w(k,j,i) ) &
[888]	1862	- ( 8.0 * ibit35 * adv_mom_5 &
	1863	+ ibit34 * adv_mom_3 &
[861]	1864	) * &
	1865	( w(k_pp,j,i) + w(k-1,j,i) ) &
[888]	1866	+ ( ibit35 * adv_mom_5 &
[861]	1867	) * &
	1868	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
	1869	)
	1870
	1871	diss_t(k) = - ABS( w_comp ) * ( &
[888]	1872	( 10.0 * ibit35 * adv_mom_5 &
	1873	+ 3.0 * ibit34 * adv_mom_3 &
	1874	+ ibit33 * adv_mom_1 &
[861]	1875	) * &
	1876	( w(k+1,j,i) - w(k,j,i) ) &
[888]	1877	- ( 5.0 * ibit35 * adv_mom_5 &
	1878	+ ibit34 * adv_mom_3 &
[861]	1879	) * &
	1880	( w(k_pp,j,i) - w(k-1,j,i) ) &
[888]	1881	+ ( ibit35 * adv_mom_5 &
[861]	1882	) * &
	1883	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
	1884	)
[862]	1885
[861]	1886	!
	1887	!-- Calculate the divergence of the velocity field. A respective
	1888	!-- correction is needed to overcome numerical instabilities introduced
	1889	!-- by a not sufficient reduction of divergences near topography.
	1890	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
	1891	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
	1892	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) &
	1893	) * 0.5
	1894
[713]	1895	tend(k,j,i) = tend(k,j,i) - ( &
	1896	( flux_r(k) + diss_r(k) &
[861]	1897	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
[713]	1898	+ ( flux_n(k) + diss_n(k) &
[861]	1899	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
	1900	+ ( flux_t(k) + diss_t(k) &
	1901	- flux_d - diss_d ) * ddzu(k+1) &
	1902	) + div * w(k,j,i)
[713]	1903
[736]	1904	flux_l_w(k,j,tn) = flux_r(k)
	1905	diss_l_w(k,j,tn) = diss_r(k)
[861]	1906	flux_s_w(k,tn) = flux_n(k)
	1907	diss_s_w(k,tn) = diss_n(k)
	1908	flux_d = flux_t(k)
	1909	diss_d = diss_t(k)
	1910	!
	1911	!-- Statistical Evaluation of w'w'.
	1912	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
	1913	+ ( flux_t(k) + diss_t(k) ) &
	1914	* weight_substep(intermediate_timestep_count)
	1915
[713]	1916	ENDDO
	1917
[861]	1918	DO k = nzb_max+1, nzt
	1919
	1920	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
	1921	flux_r(k) = u_comp * ( &
	1922	37.0 * ( w(k,j,i+1) + w(k,j,i) ) &
	1923	- 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) &
	1924	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
	1925
	1926	diss_r(k) = - ABS( u_comp ) * ( &
	1927	10.0 * ( w(k,j,i+1) - w(k,j,i) ) &
	1928	- 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) &
	1929	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
	1930
	1931	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
	1932	flux_n(k) = v_comp * ( &
	1933	37.0 * ( w(k,j+1,i) + w(k,j,i) ) &
	1934	- 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) &
	1935	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
	1936
	1937	diss_n(k) = - ABS( v_comp ) * ( &
	1938	10.0 * ( w(k,j+1,i) - w(k,j,i) ) &
	1939	- 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) &
	1940	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
[713]	1941	!
[861]	1942	!-- k index has to be modified near bottom and top, else array
	1943	!-- subscripts will be exceeded.
[1221]	1944	ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
	1945	ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
	1946	ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
[713]	1947
[888]	1948	k_ppp = k + 3 * ibit35
	1949	k_pp = k + 2 * ( 1 - ibit33 )
	1950	k_mm = k - 2 * ibit35
	1951
[713]	1952	w_comp = w(k+1,j,i) + w(k,j,i)
[861]	1953	flux_t(k) = w_comp * ( &
[888]	1954	( 37.0 * ibit35 * adv_mom_5 &
	1955	+ 7.0 * ibit34 * adv_mom_3 &
	1956	+ ibit33 * adv_mom_1 &
[861]	1957	) * &
	1958	( w(k+1,j,i) + w(k,j,i) ) &
[888]	1959	- ( 8.0 * ibit35 * adv_mom_5 &
	1960	+ ibit34 * adv_mom_3 &
[861]	1961	) * &
	1962	( w(k_pp,j,i) + w(k-1,j,i) ) &
[888]	1963	+ ( ibit35 * adv_mom_5 &
[861]	1964	) * &
	1965	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
	1966	)
[713]	1967
[861]	1968	diss_t(k) = - ABS( w_comp ) * ( &
[888]	1969	( 10.0 * ibit35 * adv_mom_5 &
	1970	+ 3.0 * ibit34 * adv_mom_3 &
	1971	+ ibit33 * adv_mom_1 &
[861]	1972	) * &
	1973	( w(k+1,j,i) - w(k,j,i) ) &
[888]	1974	- ( 5.0 * ibit35 * adv_mom_5 &
	1975	+ ibit34 * adv_mom_3 &
[861]	1976	) * &
	1977	( w(k_pp,j,i) - w(k-1,j,i) ) &
[888]	1978	+ ( ibit35 * adv_mom_5 &
[861]	1979	) * &
	1980	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
	1981	)
[713]	1982	!
[861]	1983	!-- Calculate the divergence of the velocity field. A respective
	1984	!-- correction is needed to overcome numerical instabilities introduced
	1985	!-- by a not sufficient reduction of divergences near topography.
	1986	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
	1987	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
	1988	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) * ddzu(k+1) &
	1989	) * 0.5
[713]	1990
[861]	1991	tend(k,j,i) = tend(k,j,i) - ( &
	1992	( flux_r(k) + diss_r(k) &
	1993	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
	1994	+ ( flux_n(k) + diss_n(k) &
	1995	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &
	1996	+ ( flux_t(k) + diss_t(k) &
	1997	- flux_d - diss_d ) * ddzu(k+1) &
	1998	) + div * w(k,j,i)
	1999
	2000	flux_l_w(k,j,tn) = flux_r(k)
	2001	diss_l_w(k,j,tn) = diss_r(k)
	2002	flux_s_w(k,tn) = flux_n(k)
	2003	diss_s_w(k,tn) = diss_n(k)
	2004	flux_d = flux_t(k)
	2005	diss_d = diss_t(k)
	2006	!
	2007	!-- Statistical Evaluation of w'w'.
	2008	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
	2009	+ ( flux_t(k) + diss_t(k) ) &
	2010	* weight_substep(intermediate_timestep_count)
	2011
[713]	2012	ENDDO
	2013
[861]	2014
[713]	2015	END SUBROUTINE advec_w_ws_ij
	2016
	2017
	2018	!------------------------------------------------------------------------------!
	2019	! Scalar advection - Call for all grid points
	2020	!------------------------------------------------------------------------------!
	2021	SUBROUTINE advec_s_ws( sk, sk_char )
	2022
	2023	USE arrays_3d
	2024	USE constants
	2025	USE control_parameters
	2026	USE grid_variables
	2027	USE indices
	2028	USE statistics
	2029
	2030	IMPLICIT NONE
	2031
[888]	2032	INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, &
	2033	ibit7, ibit8, j, k, k_mm, k_pp, k_ppp, tn = 0
[1010]	2034	#if defined( __nopointer )
	2035	REAL, DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
	2036	#else
[713]	2037	REAL, DIMENSION(:,:,:), POINTER :: sk
[1010]	2038	#endif
	2039	REAL :: diss_d, div, flux_d, u_comp, v_comp
	2040	REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
	2041	flux_t
	2042	REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local, swap_flux_y_local
	2043	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local, &
	2044	swap_flux_x_local
	2045	CHARACTER (LEN = *), INTENT(IN) :: sk_char
[713]	2046
	2047	!
	2048	!-- Compute the fluxes for the whole left boundary of the processor domain.
	2049	i = nxl
	2050	DO j = nys, nyn
[861]	2051
	2052	DO k = nzb+1, nzb_max
	2053
[888]	2054	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
	2055	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
	2056	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
	2057
[861]	2058	u_comp = u(k,j,i) - u_gtrans
	2059	swap_flux_x_local(k,j) = u_comp * ( &
[888]	2060	( 37.0 * ibit2 * adv_sca_5 &
	2061	+ 7.0 * ibit1 * adv_sca_3 &
	2062	+ ibit0 * adv_sca_1 &
	2063	) * &
	2064	( sk(k,j,i) + sk(k,j,i-1) ) &
	2065	- ( 8.0 * ibit2 * adv_sca_5 &
	2066	+ ibit1 * adv_sca_3 &
	2067	) * &
	2068	( sk(k,j,i+1) + sk(k,j,i-2) ) &
	2069	+ ( ibit2 * adv_sca_5 &
	2070	) * &
	2071	( sk(k,j,i+2) + sk(k,j,i-3) ) &
	2072	)
[861]	2073
	2074	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
[888]	2075	( 10.0 * ibit2 * adv_sca_5 &
	2076	+ 3.0 * ibit1 * adv_sca_3 &
	2077	+ ibit0 * adv_sca_1 &
	2078	) * &
	2079	( sk(k,j,i) - sk(k,j,i-1) ) &
	2080	- ( 5.0 * ibit2 * adv_sca_5 &
	2081	+ ibit1 * adv_sca_3 &
	2082	) * &
	2083	( sk(k,j,i+1) - sk(k,j,i-2) ) &
	2084	+ ( ibit2 * adv_sca_5 &
	2085	) * &
	2086	( sk(k,j,i+2) - sk(k,j,i-3) ) &
	2087	)
[861]	2088
	2089	ENDDO
	2090
	2091	DO k = nzb_max+1, nzt
	2092
	2093	u_comp = u(k,j,i) - u_gtrans
	2094	swap_flux_x_local(k,j) = u_comp * ( &
	2095	37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) &
	2096	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
	2097	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
[888]	2098	) * adv_sca_5
[861]	2099
	2100	swap_diss_x_local(k,j) = -ABS( u_comp ) * ( &
	2101	10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) &
	2102	- 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
	2103	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
[888]	2104	) * adv_sca_5
[861]	2105
	2106	ENDDO
	2107
[713]	2108	ENDDO
[861]	2109
[713]	2110	DO i = nxl, nxr
[861]	2111
[713]	2112	j = nys
[861]	2113	DO k = nzb+1, nzb_max
[713]	2114
[888]	2115	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
	2116	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
	2117	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
	2118
	2119	v_comp = v(k,j,i) - v_gtrans
[861]	2120	swap_flux_y_local(k) = v_comp * ( &
[888]	2121	( 37.0 * ibit5 * adv_sca_5 &
	2122	+ 7.0 * ibit4 * adv_sca_3 &
	2123	+ ibit3 * adv_sca_1 &
	2124	) * &
	2125	( sk(k,j,i) + sk(k,j-1,i) ) &
	2126	- ( 8.0 * ibit5 * adv_sca_5 &
	2127	+ ibit4 * adv_sca_3 &
	2128	) * &
	2129	( sk(k,j+1,i) + sk(k,j-2,i) ) &
	2130	+ ( ibit5 * adv_sca_5 &
	2131	) * &
	2132	( sk(k,j+2,i) + sk(k,j-3,i) ) &
[861]	2133	)
	2134
	2135	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
[888]	2136	( 10.0 * ibit5 * adv_sca_5 &
	2137	+ 3.0 * ibit4 * adv_sca_3 &
	2138	+ ibit3 * adv_sca_1 &
	2139	) * &
	2140	( sk(k,j,i) - sk(k,j-1,i) ) &
	2141	- ( 5.0 * ibit5 * adv_sca_5 &
	2142	+ ibit4 * adv_sca_3 &
	2143	) * &
	2144	( sk(k,j+1,i) - sk(k,j-2,i) ) &
	2145	+ ( ibit5 * adv_sca_5 &
	2146	) * &
	2147	( sk(k,j+2,i) - sk(k,j-3,i) ) &
[861]	2148	)
	2149
	2150	ENDDO
	2151	!
	2152	!-- Above to the top of the highest topography. No degradation necessary.
	2153	DO k = nzb_max+1, nzt
	2154
	2155	v_comp = v(k,j,i) - v_gtrans
	2156	swap_flux_y_local(k) = v_comp * ( &
	2157	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
	2158	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
	2159	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
	2160	) * adv_sca_5
	2161	swap_diss_y_local(k) = -ABS( v_comp ) * ( &
	2162	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
	2163	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
	2164	+ sk(k,j+2,i) - sk(k,j-3,i) &
	2165	) * adv_sca_5
	2166
	2167	ENDDO
	2168
[713]	2169	DO j = nys, nyn
[861]	2170
	2171	flux_t(0) = 0.0
	2172	diss_t(0) = 0.0
	2173	flux_d = 0.0
	2174	diss_d = 0.0
	2175
	2176	DO k = nzb+1, nzb_max
	2177
[888]	2178	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
	2179	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
	2180	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
	2181
[861]	2182	u_comp = u(k,j,i+1) - u_gtrans
	2183	flux_r(k) = u_comp * ( &
[888]	2184	( 37.0 * ibit2 * adv_sca_5 &
	2185	+ 7.0 * ibit1 * adv_sca_3 &
	2186	+ ibit0 * adv_sca_1 &
	2187	) * &
[861]	2188	( sk(k,j,i+1) + sk(k,j,i) ) &
[888]	2189	- ( 8.0 * ibit2 * adv_sca_5 &
	2190	+ ibit1 * adv_sca_3 &
	2191	) * &
[861]	2192	( sk(k,j,i+2) + sk(k,j,i-1) ) &
[888]	2193	+ ( ibit2 * adv_sca_5 &
	2194	) * &
	2195	( sk(k,j,i+3) + sk(k,j,i-2) ) &
[861]	2196	)
	2197
	2198	diss_r(k) = -ABS( u_comp ) * ( &
[888]	2199	( 10.0 * ibit2 * adv_sca_5 &
	2200	+ 3.0 * ibit1 * adv_sca_3 &
	2201	+ ibit0 * adv_sca_1 &
	2202	) * &
[861]	2203	( sk(k,j,i+1) - sk(k,j,i) ) &
[888]	2204	- ( 5.0 * ibit2 * adv_sca_5 &
	2205	+ ibit1 * adv_sca_3 &
	2206	) * &
[861]	2207	( sk(k,j,i+2) - sk(k,j,i-1) ) &
[888]	2208	+ ( ibit2 * adv_sca_5 &
	2209	) * &
[861]	2210	( sk(k,j,i+3) - sk(k,j,i-2) ) &
	2211	)
	2212
[888]	2213	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
	2214	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
	2215	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
	2216
[861]	2217	v_comp = v(k,j+1,i) - v_gtrans
	2218	flux_n(k) = v_comp * ( &
[888]	2219	( 37.0 * ibit5 * adv_sca_5 &
	2220	+ 7.0 * ibit4 * adv_sca_3 &
	2221	+ ibit3 * adv_sca_1 &
	2222	) * &
[861]	2223	( sk(k,j+1,i) + sk(k,j,i) ) &
[888]	2224	- ( 8.0 * ibit5 * adv_sca_5 &
	2225	+ ibit4 * adv_sca_3 &
	2226	) * &
[861]	2227	( sk(k,j+2,i) + sk(k,j-1,i) ) &
[888]	2228	+ ( ibit5 * adv_sca_5 &
	2229	) * &
[861]	2230	( sk(k,j+3,i) + sk(k,j-2,i) ) &
	2231	)
	2232
	2233	diss_n(k) = -ABS( v_comp ) * ( &
[888]	2234	( 10.0 * ibit5 * adv_sca_5 &
	2235	+ 3.0 * ibit4 * adv_sca_3 &
	2236	+ ibit3 * adv_sca_1 &
	2237	) * &
[861]	2238	( sk(k,j+1,i) - sk(k,j,i) ) &
[888]	2239	- ( 5.0 * ibit5 * adv_sca_5 &
	2240	+ ibit4 * adv_sca_3 &
	2241	) * &
[861]	2242	( sk(k,j+2,i) - sk(k,j-1,i) ) &
[888]	2243	+ ( ibit5 * adv_sca_5 &
	2244	) * &
[861]	2245	( sk(k,j+3,i) - sk(k,j-2,i) ) &
	2246	)
[713]	2247	!
[861]	2248	!-- k index has to be modified near bottom and top, else array
	2249	!-- subscripts will be exceeded.
[888]	2250	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
	2251	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
	2252	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
[861]	2253
[888]	2254	k_ppp = k + 3 * ibit8
	2255	k_pp = k + 2 * ( 1 - ibit6 )
	2256	k_mm = k - 2 * ibit8
[861]	2257
[888]	2258
	2259	flux_t(k) = w(k,j,i) * ( &
	2260	( 37.0 * ibit8 * adv_sca_5 &
	2261	+ 7.0 * ibit7 * adv_sca_3 &
	2262	+ ibit6 * adv_sca_1 &
	2263	) * &
	2264	( sk(k+1,j,i) + sk(k,j,i) ) &
	2265	- ( 8.0 * ibit8 * adv_sca_5 &
	2266	+ ibit7 * adv_sca_3 &
	2267	) * &
	2268	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
	2269	+ ( ibit8 * adv_sca_5 &
	2270	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
[861]	2271	)
	2272
[888]	2273	diss_t(k) = -ABS( w(k,j,i) ) * ( &
	2274	( 10.0 * ibit8 * adv_sca_5 &
	2275	+ 3.0 * ibit7 * adv_sca_3 &
	2276	+ ibit6 * adv_sca_1 &
	2277	) * &
	2278	( sk(k+1,j,i) - sk(k,j,i) ) &
	2279	- ( 5.0 * ibit8 * adv_sca_5 &
	2280	+ ibit7 * adv_sca_3 &
	2281	) * &
	2282	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
	2283	+ ( ibit8 * adv_sca_5 &
	2284	) * &
	2285	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
	2286	)
[861]	2287	!
	2288	!-- Calculate the divergence of the velocity field. A respective
	2289	!-- correction is needed to overcome numerical instabilities caused
	2290	!-- by a not sufficient reduction of divergences near topography.
	2291	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
	2292	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
	2293	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	2294
	2295	tend(k,j,i) = tend(k,j,i) - ( &
	2296	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
	2297	swap_diss_x_local(k,j) ) * ddx &
	2298	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
	2299	swap_diss_y_local(k) ) * ddy &
	2300	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
	2301	) * ddzw(k) &
	2302	) + sk(k,j,i) * div
	2303
	2304	swap_flux_y_local(k) = flux_n(k)
	2305	swap_diss_y_local(k) = diss_n(k)
[713]	2306	swap_flux_x_local(k,j) = flux_r(k)
	2307	swap_diss_x_local(k,j) = diss_r(k)
[861]	2308	flux_d = flux_t(k)
	2309	diss_d = diss_t(k)
[713]	2310
[861]	2311	ENDDO
	2312
	2313	DO k = nzb_max+1, nzt
	2314
	2315	u_comp = u(k,j,i+1) - u_gtrans
	2316	flux_r(k) = u_comp * ( &
	2317	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
	2318	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
	2319	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
	2320	diss_r(k) = -ABS( u_comp ) * ( &
	2321	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
	2322	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
	2323	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
	2324
	2325	v_comp = v(k,j+1,i) - v_gtrans
	2326	flux_n(k) = v_comp * ( &
	2327	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
	2328	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
	2329	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
	2330	diss_n(k) = -ABS( v_comp ) * ( &
	2331	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
	2332	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
	2333	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
[713]	2334	!
[861]	2335	!-- k index has to be modified near bottom and top, else array
	2336	!-- subscripts will be exceeded.
[888]	2337	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
	2338	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
	2339	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
[713]	2340
[888]	2341	k_ppp = k + 3 * ibit8
	2342	k_pp = k + 2 * ( 1 - ibit6 )
	2343	k_mm = k - 2 * ibit8
	2344
	2345
	2346	flux_t(k) = w(k,j,i) * ( &
	2347	( 37.0 * ibit8 * adv_sca_5 &
	2348	+ 7.0 * ibit7 * adv_sca_3 &
	2349	+ ibit6 * adv_sca_1 &
	2350	) * &
	2351	( sk(k+1,j,i) + sk(k,j,i) ) &
	2352	- ( 8.0 * ibit8 * adv_sca_5 &
	2353	+ ibit7 * adv_sca_3 &
	2354	) * &
	2355	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
	2356	+ ( ibit8 * adv_sca_5 &
	2357	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
[861]	2358	)
	2359
[888]	2360	diss_t(k) = -ABS( w(k,j,i) ) * ( &
	2361	( 10.0 * ibit8 * adv_sca_5 &
	2362	+ 3.0 * ibit7 * adv_sca_3 &
	2363	+ ibit6 * adv_sca_1 &
	2364	) * &
	2365	( sk(k+1,j,i) - sk(k,j,i) ) &
	2366	- ( 5.0 * ibit8 * adv_sca_5 &
	2367	+ ibit7 * adv_sca_3 &
	2368	) * &
	2369	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
	2370	+ ( ibit8 * adv_sca_5 &
	2371	) * &
	2372	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
	2373	)
[713]	2374	!
[861]	2375	!-- Calculate the divergence of the velocity field. A respective
	2376	!-- correction is needed to overcome numerical instabilities introduced
	2377	!-- by a not sufficient reduction of divergences near topography.
	2378	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
	2379	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
	2380	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
[713]	2381
[861]	2382	tend(k,j,i) = tend(k,j,i) - ( &
	2383	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j) - &
	2384	swap_diss_x_local(k,j) ) * ddx &
	2385	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k) - &
	2386	swap_diss_y_local(k) ) * ddy &
	2387	+ ( flux_t(k) + diss_t(k) - flux_d - diss_d &
	2388	) * ddzw(k) &
	2389	) + sk(k,j,i) * div
	2390
	2391	swap_flux_y_local(k) = flux_n(k)
	2392	swap_diss_y_local(k) = diss_n(k)
	2393	swap_flux_x_local(k,j) = flux_r(k)
	2394	swap_diss_x_local(k,j) = diss_r(k)
	2395	flux_d = flux_t(k)
	2396	diss_d = diss_t(k)
	2397
	2398	ENDDO
[713]	2399	!
[861]	2400	!-- evaluation of statistics
	2401	SELECT CASE ( sk_char )
[713]	2402
[861]	2403	CASE ( 'pt' )
	2404	DO k = nzb, nzt
	2405	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
	2406	+ ( flux_t(k) + diss_t(k) ) &
	2407	* weight_substep(intermediate_timestep_count)
	2408	ENDDO
	2409	CASE ( 'sa' )
	2410	DO k = nzb, nzt
	2411	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
	2412	+ ( flux_t(k) + diss_t(k) ) &
	2413	* weight_substep(intermediate_timestep_count)
	2414	ENDDO
	2415	CASE ( 'q' )
	2416	DO k = nzb, nzt
	2417	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
	2418	+ ( flux_t(k) + diss_t(k) ) &
	2419	* weight_substep(intermediate_timestep_count)
	2420	ENDDO
[713]	2421
[861]	2422	END SELECT
	2423
[713]	2424	ENDDO
	2425	ENDDO
	2426
	2427	END SUBROUTINE advec_s_ws
	2428
	2429
	2430	!------------------------------------------------------------------------------!
[1015]	2431	! Scalar advection - Call for all grid points - accelerator version
	2432	!------------------------------------------------------------------------------!
	2433	SUBROUTINE advec_s_ws_acc ( sk, sk_char )
	2434
	2435	USE arrays_3d
	2436	USE constants
	2437	USE control_parameters
	2438	USE grid_variables
	2439	USE indices
	2440	USE statistics
	2441
	2442	IMPLICIT NONE
	2443
	2444	CHARACTER (LEN = *), INTENT(IN) :: sk_char
	2445
	2446	INTEGER :: i, ibit0, ibit1, ibit2, ibit3, ibit4, ibit5, ibit6, &
	2447	ibit7, ibit8, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
	2448
	2449	REAL :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, flux_d, &
	2450	flux_l, flux_n, flux_r, flux_s, flux_t, u_comp, v_comp
	2451
	2452	REAL, INTENT(IN), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk
	2453
	2454	!
	2455	!-- Computation of fluxes and tendency terms
[1221]	2456	!$acc kernels present( ddzw, sk, tend, u, v, w, wall_flags_0, wall_flags_00 )
[1015]	2457	!$acc loop
[1128]	2458	DO i = i_left, i_right
	2459	DO j = j_south, j_north
[1015]	2460	!$acc loop vector( 32 )
	2461	DO k = nzb+1, nzt
	2462
	2463	ibit2 = IBITS(wall_flags_0(k,j,i),2,1)
	2464	ibit1 = IBITS(wall_flags_0(k,j,i),1,1)
	2465	ibit0 = IBITS(wall_flags_0(k,j,i),0,1)
	2466
	2467	u_comp = u(k,j,i) - u_gtrans
	2468	flux_l = u_comp * ( &
	2469	( 37.0 * ibit2 * adv_sca_5 &
	2470	+ 7.0 * ibit1 * adv_sca_3 &
	2471	+ ibit0 * adv_sca_1 &
	2472	) * &
	2473	( sk(k,j,i) + sk(k,j,i-1) ) &
	2474	- ( 8.0 * ibit2 * adv_sca_5 &
	2475	+ ibit1 * adv_sca_3 &
	2476	) * &
	2477	( sk(k,j,i+1) + sk(k,j,i-2) ) &
	2478	+ ( ibit2 * adv_sca_5 &
	2479	) * &
	2480	( sk(k,j,i+2) + sk(k,j,i-3) ) &
	2481	)
	2482
	2483	diss_l = -ABS( u_comp ) * ( &
	2484	( 10.0 * ibit2 * adv_sca_5 &
	2485	+ 3.0 * ibit1 * adv_sca_3 &
	2486	+ ibit0 * adv_sca_1 &
	2487	) * &
	2488	( sk(k,j,i) - sk(k,j,i-1) ) &
	2489	- ( 5.0 * ibit2 * adv_sca_5 &
	2490	+ ibit1 * adv_sca_3 &
	2491	) * &
	2492	( sk(k,j,i+1) - sk(k,j,i-2) ) &
	2493	+ ( ibit2 * adv_sca_5 &
	2494	) * &
	2495	( sk(k,j,i+2) - sk(k,j,i-3) ) &
	2496	)
	2497
	2498	u_comp = u(k,j,i+1) - u_gtrans
	2499	flux_r = u_comp * ( &
	2500	( 37.0 * ibit2 * adv_sca_5 &
	2501	+ 7.0 * ibit1 * adv_sca_3 &
	2502	+ ibit0 * adv_sca_1 &
	2503	) * &
	2504	( sk(k,j,i+1) + sk(k,j,i) ) &
	2505	- ( 8.0 * ibit2 * adv_sca_5 &
	2506	+ ibit1 * adv_sca_3 &
	2507	) * &
	2508	( sk(k,j,i+2) + sk(k,j,i-1) ) &
	2509	+ ( ibit2 * adv_sca_5 &
	2510	) * &
	2511	( sk(k,j,i+3) + sk(k,j,i-2) ) &
	2512	)
	2513
	2514	diss_r = -ABS( u_comp ) * ( &
	2515	( 10.0 * ibit2 * adv_sca_5 &
	2516	+ 3.0 * ibit1 * adv_sca_3 &
	2517	+ ibit0 * adv_sca_1 &
	2518	) * &
	2519	( sk(k,j,i+1) - sk(k,j,i) ) &
	2520	- ( 5.0 * ibit2 * adv_sca_5 &
	2521	+ ibit1 * adv_sca_3 &
	2522	) * &
	2523	( sk(k,j,i+2) - sk(k,j,i-1) ) &
	2524	+ ( ibit2 * adv_sca_5 &
	2525	) * &
	2526	( sk(k,j,i+3) - sk(k,j,i-2) ) &
	2527	)
	2528
	2529	ibit5 = IBITS(wall_flags_0(k,j,i),5,1)
	2530	ibit4 = IBITS(wall_flags_0(k,j,i),4,1)
	2531	ibit3 = IBITS(wall_flags_0(k,j,i),3,1)
	2532
	2533	v_comp = v(k,j,i) - v_gtrans
	2534	flux_s = v_comp * ( &
	2535	( 37.0 * ibit5 * adv_sca_5 &
	2536	+ 7.0 * ibit4 * adv_sca_3 &
	2537	+ ibit3 * adv_sca_1 &
	2538	) * &
	2539	( sk(k,j,i) + sk(k,j-1,i) ) &
	2540	- ( 8.0 * ibit5 * adv_sca_5 &
	2541	+ ibit4 * adv_sca_3 &
	2542	) * &
	2543	( sk(k,j+1,i) + sk(k,j-2,i) ) &
	2544	+ ( ibit5 * adv_sca_5 &
	2545	) * &
	2546	( sk(k,j+2,i) + sk(k,j-3,i) ) &
	2547	)
	2548
	2549	diss_s = -ABS( v_comp ) * ( &
	2550	( 10.0 * ibit5 * adv_sca_5 &
	2551	+ 3.0 * ibit4 * adv_sca_3 &
	2552	+ ibit3 * adv_sca_1 &
	2553	) * &
	2554	( sk(k,j,i) - sk(k,j-1,i) ) &
	2555	- ( 5.0 * ibit5 * adv_sca_5 &
	2556	+ ibit4 * adv_sca_3 &
	2557	) * &
	2558	( sk(k,j+1,i) - sk(k,j-2,i) ) &
	2559	+ ( ibit5 * adv_sca_5 &
	2560	) * &
	2561	( sk(k,j+2,i) - sk(k,j-3,i) ) &
	2562	)
	2563
	2564
	2565	v_comp = v(k,j+1,i) - v_gtrans
	2566	flux_n = v_comp * ( &
	2567	( 37.0 * ibit5 * adv_sca_5 &
	2568	+ 7.0 * ibit4 * adv_sca_3 &
	2569	+ ibit3 * adv_sca_1 &
	2570	) * &
	2571	( sk(k,j+1,i) + sk(k,j,i) ) &
	2572	- ( 8.0 * ibit5 * adv_sca_5 &
	2573	+ ibit4 * adv_sca_3 &
	2574	) * &
	2575	( sk(k,j+2,i) + sk(k,j-1,i) ) &
	2576	+ ( ibit5 * adv_sca_5 &
	2577	) * &
	2578	( sk(k,j+3,i) + sk(k,j-2,i) ) &
	2579	)
	2580
	2581	diss_n = -ABS( v_comp ) * ( &
	2582	( 10.0 * ibit5 * adv_sca_5 &
	2583	+ 3.0 * ibit4 * adv_sca_3 &
	2584	+ ibit3 * adv_sca_1 &
	2585	) * &
	2586	( sk(k,j+1,i) - sk(k,j,i) ) &
	2587	- ( 5.0 * ibit5 * adv_sca_5 &
	2588	+ ibit4 * adv_sca_3 &
	2589	) * &
	2590	( sk(k,j+2,i) - sk(k,j-1,i) ) &
	2591	+ ( ibit5 * adv_sca_5 &
	2592	) * &
	2593	( sk(k,j+3,i) - sk(k,j-2,i) ) &
	2594	)
	2595
	2596	!
	2597	!-- indizes k_m, k_mm, ... should be known at these point
	2598	ibit8 = IBITS(wall_flags_0(k-1,j,i),8,1)
	2599	ibit7 = IBITS(wall_flags_0(k-1,j,i),7,1)
	2600	ibit6 = IBITS(wall_flags_0(k-1,j,i),6,1)
	2601
[1027]	2602	k_pp = k + 2 * ibit8
	2603	k_mm = k - 2 * ( ibit7 + ibit8 )
[1015]	2604	k_mmm = k - 3 * ibit8
	2605
	2606	flux_d = w(k-1,j,i) * ( &
	2607	( 37.0 * ibit8 * adv_sca_5 &
	2608	+ 7.0 * ibit7 * adv_sca_3 &
	2609	+ ibit6 * adv_sca_1 &
	2610	) * &
	2611	( sk(k,j,i) + sk(k-1,j,i) ) &
	2612	- ( 8.0 * ibit8 * adv_sca_5 &
	2613	+ ibit7 * adv_sca_3 &
	2614	) * &
	2615	( sk(k+1,j,i) + sk(k_mm,j,i) ) &
	2616	+ ( ibit8 * adv_sca_5 &
	2617	) * ( sk(k_pp,j,i)+ sk(k_mmm,j,i) ) &
	2618	)
	2619
	2620	diss_d = -ABS( w(k-1,j,i) ) * ( &
	2621	( 10.0 * ibit8 * adv_sca_5 &
	2622	+ 3.0 * ibit7 * adv_sca_3 &
	2623	+ ibit6 * adv_sca_1 &
	2624	) * &
	2625	( sk(k,j,i) - sk(k-1,j,i) ) &
	2626	- ( 5.0 * ibit8 * adv_sca_5 &
	2627	+ ibit7 * adv_sca_3 &
	2628	) * &
	2629	( sk(k+1,j,i) - sk(k_mm,j,i) ) &
	2630	+ ( ibit8 * adv_sca_5 &
	2631	) * &
	2632	( sk(k_pp,j,i) - sk(k_mmm,j,i) ) &
	2633	)
	2634
	2635	ibit8 = IBITS(wall_flags_0(k,j,i),8,1)
	2636	ibit7 = IBITS(wall_flags_0(k,j,i),7,1)
	2637	ibit6 = IBITS(wall_flags_0(k,j,i),6,1)
	2638
	2639	k_ppp = k + 3 * ibit8
	2640	k_pp = k + 2 * ( 1 - ibit6 )
	2641	k_mm = k - 2 * ibit8
	2642
	2643	flux_t = w(k,j,i) * ( &
	2644	( 37.0 * ibit8 * adv_sca_5 &
	2645	+ 7.0 * ibit7 * adv_sca_3 &
	2646	+ ibit6 * adv_sca_1 &
	2647	) * &
	2648	( sk(k+1,j,i) + sk(k,j,i) ) &
	2649	- ( 8.0 * ibit8 * adv_sca_5 &
	2650	+ ibit7 * adv_sca_3 &
	2651	) * &
	2652	( sk(k_pp,j,i) + sk(k-1,j,i) ) &
	2653	+ ( ibit8 * adv_sca_5 &
	2654	) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &
	2655	)
	2656
	2657	diss_t = -ABS( w(k,j,i) ) * ( &
	2658	( 10.0 * ibit8 * adv_sca_5 &
	2659	+ 3.0 * ibit7 * adv_sca_3 &
	2660	+ ibit6 * adv_sca_1 &
	2661	) * &
	2662	( sk(k+1,j,i) - sk(k,j,i) ) &
	2663	- ( 5.0 * ibit8 * adv_sca_5 &
	2664	+ ibit7 * adv_sca_3 &
	2665	) * &
	2666	( sk(k_pp,j,i) - sk(k-1,j,i) ) &
	2667	+ ( ibit8 * adv_sca_5 &
	2668	) * &
	2669	( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &
	2670	)
	2671	!
	2672	!-- Calculate the divergence of the velocity field. A respective
	2673	!-- correction is needed to overcome numerical instabilities caused
	2674	!-- by a not sufficient reduction of divergences near topography.
	2675	div = ( u(k,j,i+1) - u(k,j,i) ) * ddx &
	2676	+ ( v(k,j+1,i) - v(k,j,i) ) * ddy &
	2677	+ ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
	2678
	2679	tend(k,j,i) = - ( &
	2680	( flux_r + diss_r - flux_l - diss_l ) * ddx &
	2681	+ ( flux_n + diss_n - flux_s - diss_s ) * ddy &
	2682	+ ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k)&
	2683	) + div * sk(k,j,i)
	2684
	2685	!++
	2686	!-- Evaluation of statistics
	2687	! SELECT CASE ( sk_char )
	2688	!
	2689	! CASE ( 'pt' )
	2690	! sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
	2691	! + ( flux_t + diss_t ) &
	2692	! * weight_substep(intermediate_timestep_count)
	2693	! CASE ( 'sa' )
	2694	! sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
	2695	! + ( flux_t + diss_t ) &
	2696	! * weight_substep(intermediate_timestep_count)
	2697	! CASE ( 'q' )
	2698	! sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
	2699	! + ( flux_t + diss_t ) &
	2700	! * weight_substep(intermediate_timestep_count)
	2701	!
	2702	! END SELECT
	2703
	2704	ENDDO
	2705	ENDDO
	2706	ENDDO
	2707	!$acc end kernels
	2708
	2709	END SUBROUTINE advec_s_ws_acc
	2710
	2711
	2712	!------------------------------------------------------------------------------!
[713]	2713	! Advection of u - Call for all grid points
	2714	!------------------------------------------------------------------------------!
	2715	SUBROUTINE advec_u_ws
	2716
	2717	USE arrays_3d
	2718	USE constants
	2719	USE control_parameters
	2720	USE grid_variables
	2721	USE indices
	2722	USE statistics
	2723
	2724	IMPLICIT NONE
	2725
[888]	2726	INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, &
	2727	ibit16, ibit17, j, k, k_mm, k_pp, k_ppp, tn = 0
[862]	2728	REAL :: diss_d, div, flux_d, gu, gv, v_comp, w_comp
	2729	REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_u, swap_flux_y_local_u
	2730	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_u, &
	2731	swap_flux_x_local_u
	2732	REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
	2733	flux_t, u_comp
[713]	2734
	2735	gu = 2.0 * u_gtrans
	2736	gv = 2.0 * v_gtrans
	2737
	2738	!
	2739	!-- Compute the fluxes for the whole left boundary of the processor domain.
	2740	i = nxlu
	2741	DO j = nys, nyn
[861]	2742	DO k = nzb+1, nzb_max
	2743
[888]	2744	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
	2745	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
	2746	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
	2747
[861]	2748	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
	2749	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
[888]	2750	( 37.0 * ibit11 * adv_mom_5 &
	2751	+ 7.0 * ibit10 * adv_mom_3 &
	2752	+ ibit9 * adv_mom_1 &
	2753	) * &
[861]	2754	( u(k,j,i) + u(k,j,i-1) ) &
[888]	2755	- ( 8.0 * ibit11 * adv_mom_5 &
	2756	+ ibit10 * adv_mom_3 &
	2757	) * &
[861]	2758	( u(k,j,i+1) + u(k,j,i-2) ) &
[888]	2759	+ ( ibit11 * adv_mom_5 &
	2760	) * &
[861]	2761	( u(k,j,i+2) + u(k,j,i-3) ) &
	2762	)
	2763
	2764	swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( &
[888]	2765	( 10.0 * ibit11 * adv_mom_5 &
	2766	+ 3.0 * ibit10 * adv_mom_3 &
	2767	+ ibit9 * adv_mom_1 &
	2768	) * &
[861]	2769	( u(k,j,i) - u(k,j,i-1) ) &
[888]	2770	- ( 5.0 * ibit11 * adv_mom_5 &
	2771	+ ibit10 * adv_mom_3 &
	2772	) * &
[861]	2773	( u(k,j,i+1) - u(k,j,i-2) ) &
[888]	2774	+ ( ibit11 * adv_mom_5 &
	2775	) * &
[861]	2776	( u(k,j,i+2) - u(k,j,i-3) ) &
	2777	)
	2778
	2779	ENDDO
	2780
	2781	DO k = nzb_max+1, nzt
	2782
	2783	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
	2784	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
	2785	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
	2786	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
	2787	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
	2788	swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( &
	2789	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
	2790	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
	2791	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
	2792
	2793	ENDDO
[713]	2794	ENDDO
	2795
	2796	DO i = nxlu, nxr
	2797	!
[861]	2798	!-- The following loop computes the fluxes for the south boundary points
	2799	j = nys
	2800	DO k = nzb+1, nzb_max
[713]	2801
[888]	2802	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
	2803	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
	2804	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
	2805
[861]	2806	v_comp = v(k,j,i) + v(k,j,i-1) - gv
	2807	swap_flux_y_local_u(k) = v_comp * ( &
[888]	2808	( 37.0 * ibit14 * adv_mom_5 &
	2809	+ 7.0 * ibit13 * adv_mom_3 &
	2810	+ ibit12 * adv_mom_1 &
	2811	) * &
[861]	2812	( u(k,j,i) + u(k,j-1,i) ) &
[888]	2813	- ( 8.0 * ibit14 * adv_mom_5 &
	2814	+ ibit13 * adv_mom_3 &
	2815	) * &
[861]	2816	( u(k,j+1,i) + u(k,j-2,i) ) &
[888]	2817	+ ( ibit14 * adv_mom_5 &
	2818	) * &
[861]	2819	( u(k,j+2,i) + u(k,j-3,i) ) &
	2820	)
[713]	2821
[861]	2822	swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( &
[888]	2823	( 10.0 * ibit14 * adv_mom_5 &
	2824	+ 3.0 * ibit13 * adv_mom_3 &
	2825	+ ibit12 * adv_mom_1 &
	2826	) * &
[861]	2827	( u(k,j,i) - u(k,j-1,i) ) &
[888]	2828	- ( 5.0 * ibit14 * adv_mom_5 &
	2829	+ ibit13 * adv_mom_3 &
	2830	) * &
[861]	2831	( u(k,j+1,i) - u(k,j-2,i) ) &
[888]	2832	+ ( ibit14 * adv_mom_5 &
	2833	) * &
[861]	2834	( u(k,j+2,i) - u(k,j-3,i) ) &
	2835	)
	2836
	2837	ENDDO
	2838
	2839	DO k = nzb_max+1, nzt
	2840
	2841	v_comp = v(k,j,i) + v(k,j,i-1) - gv
	2842	swap_flux_y_local_u(k) = v_comp * ( &
	2843	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
	2844	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
	2845	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
	2846	swap_diss_y_local_u(k) = - ABS(v_comp) * ( &
	2847	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
	2848	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
	2849	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
	2850
	2851	ENDDO
[713]	2852	!
[861]	2853	!-- Computation of interior fluxes and tendency terms
[713]	2854	DO j = nys, nyn
[861]	2855
	2856	flux_t(0) = 0.0
	2857	diss_t(0) = 0.0
	2858	flux_d = 0.0
	2859	diss_d = 0.0
	2860
	2861	DO k = nzb+1, nzb_max
	2862
[888]	2863	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
	2864	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
	2865	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
	2866
[861]	2867	u_comp(k) = u(k,j,i+1) + u(k,j,i)
	2868	flux_r(k) = ( u_comp(k) - gu ) * ( &
[888]	2869	( 37.0 * ibit11 * adv_mom_5 &
	2870	+ 7.0 * ibit10 * adv_mom_3 &
	2871	+ ibit9 * adv_mom_1 &
	2872	) * &
[861]	2873	( u(k,j,i+1) + u(k,j,i) ) &
[888]	2874	- ( 8.0 * ibit11 * adv_mom_5 &
	2875	+ ibit10 * adv_mom_3 &
	2876	) * &
[861]	2877	( u(k,j,i+2) + u(k,j,i-1) ) &
[888]	2878	+ ( ibit11 * adv_mom_5 &
	2879	) * &
[861]	2880	( u(k,j,i+3) + u(k,j,i-2) ) &
	2881	)
	2882
	2883	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
[888]	2884	( 10.0 * ibit11 * adv_mom_5 &
	2885	+ 3.0 * ibit10 * adv_mom_3 &
	2886	+ ibit9 * adv_mom_1 &
	2887	) * &
[861]	2888	( u(k,j,i+1) - u(k,j,i) ) &
[888]	2889	- ( 5.0 * ibit11 * adv_mom_5 &
	2890	+ ibit10 * adv_mom_3 &
	2891	) * &
[861]	2892	( u(k,j,i+2) - u(k,j,i-1) ) &
[888]	2893	+ ( ibit11 * adv_mom_5 &
	2894	) * &
[861]	2895	( u(k,j,i+3) - u(k,j,i-2) ) &
	2896	)
	2897
[888]	2898	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
	2899	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
	2900	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
	2901
[861]	2902	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
	2903	flux_n(k) = v_comp * ( &
[888]	2904	( 37.0 * ibit14 * adv_mom_5 &
	2905	+ 7.0 * ibit13 * adv_mom_3 &
	2906	+ ibit12 * adv_mom_1 &
	2907	) * &
[861]	2908	( u(k,j+1,i) + u(k,j,i) ) &
[888]	2909	- ( 8.0 * ibit14 * adv_mom_5 &
	2910	+ ibit13 * adv_mom_3 &
	2911	) * &
[861]	2912	( u(k,j+2,i) + u(k,j-1,i) ) &
[888]	2913	+ ( ibit14 * adv_mom_5 &
	2914	) * &
[861]	2915	( u(k,j+3,i) + u(k,j-2,i) ) &
	2916	)
	2917
	2918	diss_n(k) = - ABS ( v_comp ) * ( &
[888]	2919	( 10.0 * ibit14 * adv_mom_5 &
	2920	+ 3.0 * ibit13 * adv_mom_3 &
	2921	+ ibit12 * adv_mom_1 &
	2922	) * &
[861]	2923	( u(k,j+1,i) - u(k,j,i) ) &
[888]	2924	- ( 5.0 * ibit14 * adv_mom_5 &
	2925	+ ibit13 * adv_mom_3 &
	2926	) * &
[861]	2927	( u(k,j+2,i) - u(k,j-1,i) ) &
[888]	2928	+ ( ibit14 * adv_mom_5 &
	2929	) * &
[861]	2930	( u(k,j+3,i) - u(k,j-2,i) ) &
	2931	)
[713]	2932	!
[861]	2933	!-- k index has to be modified near bottom and top, else array
	2934	!-- subscripts will be exceeded.
[888]	2935	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
	2936	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
	2937	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
[861]	2938
[888]	2939	k_ppp = k + 3 * ibit17
	2940	k_pp = k + 2 * ( 1 - ibit15 )
	2941	k_mm = k - 2 * ibit17
	2942
[861]	2943	w_comp = w(k,j,i) + w(k,j,i-1)
	2944	flux_t(k) = w_comp * ( &
[888]	2945	( 37.0 * ibit17 * adv_mom_5 &
	2946	+ 7.0 * ibit16 * adv_mom_3 &
	2947	+ ibit15 * adv_mom_1 &
	2948	) * &
[861]	2949	( u(k+1,j,i) + u(k,j,i) ) &
[888]	2950	- ( 8.0 * ibit17 * adv_mom_5 &
	2951	+ ibit16 * adv_mom_3 &
	2952	) * &
[861]	2953	( u(k_pp,j,i) + u(k-1,j,i) ) &
[888]	2954	+ ( ibit17 * adv_mom_5 &
	2955	) * &
[861]	2956	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
	2957	)
	2958
	2959	diss_t(k) = - ABS( w_comp ) * ( &
[888]	2960	( 10.0 * ibit17 * adv_mom_5 &
	2961	+ 3.0 * ibit16 * adv_mom_3 &
	2962	+ ibit15 * adv_mom_1 &
	2963	) * &
[861]	2964	( u(k+1,j,i) - u(k,j,i) ) &
[888]	2965	- ( 5.0 * ibit17 * adv_mom_5 &
	2966	+ ibit16 * adv_mom_3 &
	2967	) * &
[861]	2968	( u(k_pp,j,i) - u(k-1,j,i) ) &
[888]	2969	+ ( ibit17 * adv_mom_5 &
	2970	) * &
[861]	2971	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
	2972	)
[713]	2973	!
[861]	2974	!-- Calculate the divergence of the velocity field. A respective
	2975	!-- correction is needed to overcome numerical instabilities caused
	2976	!-- by a not sufficient reduction of divergences near topography.
	2977	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
	2978	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
	2979	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) &
	2980	* ddzw(k) &
	2981	) * 0.5
[713]	2982
[861]	2983	tend(k,j,i) = tend(k,j,i) - ( &
	2984	( flux_r(k) + diss_r(k) &
	2985	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
	2986	+ ( flux_n(k) + diss_n(k) &
	2987	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
	2988	+ ( flux_t(k) + diss_t(k) &
	2989	- flux_d - diss_d &
	2990	) * ddzw(k) &
	2991	) + div * u(k,j,i)
	2992
	2993	swap_flux_x_local_u(k,j) = flux_r(k)
	2994	swap_diss_x_local_u(k,j) = diss_r(k)
	2995	swap_flux_y_local_u(k) = flux_n(k)
	2996	swap_diss_y_local_u(k) = diss_n(k)
	2997	flux_d = flux_t(k)
	2998	diss_d = diss_t(k)
[713]	2999	!
[861]	3000	!-- Statistical Evaluation of u'u'. The factor has to be applied
	3001	!-- for right evaluation when gallilei_trans = .T. .
	3002	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
	3003	+ ( flux_r(k) * &
	3004	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	3005	/ ( u_comp(k) - gu + 1.0E-20 ) &
	3006	+ diss_r(k) * &
	3007	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	3008	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
	3009	* weight_substep(intermediate_timestep_count)
	3010	!
	3011	!-- Statistical Evaluation of w'u'.
	3012	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
	3013	+ ( flux_t(k) + diss_t(k) ) &
	3014	* weight_substep(intermediate_timestep_count)
	3015	ENDDO
[713]	3016
[861]	3017	DO k = nzb_max+1, nzt
[713]	3018
[861]	3019	u_comp(k) = u(k,j,i+1) + u(k,j,i)
	3020	flux_r(k) = ( u_comp(k) - gu ) * ( &
	3021	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
	3022	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
	3023	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
	3024	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
	3025	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
	3026	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
	3027	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
[713]	3028
[861]	3029	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
	3030	flux_n(k) = v_comp * ( &
	3031	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
	3032	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
	3033	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
	3034	diss_n(k) = - ABS( v_comp ) * ( &
	3035	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
	3036	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
	3037	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
[713]	3038	!
[861]	3039	!-- k index has to be modified near bottom and top, else array
	3040	!-- subscripts will be exceeded.
[888]	3041	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
	3042	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
	3043	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
[861]	3044
[888]	3045	k_ppp = k + 3 * ibit17
	3046	k_pp = k + 2 * ( 1 - ibit15 )
	3047	k_mm = k - 2 * ibit17
	3048
[861]	3049	w_comp = w(k,j,i) + w(k,j,i-1)
	3050	flux_t(k) = w_comp * ( &
[888]	3051	( 37.0 * ibit17 * adv_mom_5 &
	3052	+ 7.0 * ibit16 * adv_mom_3 &
	3053	+ ibit15 * adv_mom_1 &
	3054	) * &
[861]	3055	( u(k+1,j,i) + u(k,j,i) ) &
[888]	3056	- ( 8.0 * ibit17 * adv_mom_5 &
	3057	+ ibit16 * adv_mom_3 &
	3058	) * &
[861]	3059	( u(k_pp,j,i) + u(k-1,j,i) ) &
[888]	3060	+ ( ibit17 * adv_mom_5 &
	3061	) * &
[861]	3062	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
	3063	)
	3064
[888]	3065	diss_t(k) = - ABS( w_comp ) * ( &
	3066	( 10.0 * ibit17 * adv_mom_5 &
	3067	+ 3.0 * ibit16 * adv_mom_3 &
	3068	+ ibit15 * adv_mom_1 &
	3069	) * &
[861]	3070	( u(k+1,j,i) - u(k,j,i) ) &
[888]	3071	- ( 5.0 * ibit17 * adv_mom_5 &
	3072	+ ibit16 * adv_mom_3 &
	3073	) * &
[861]	3074	( u(k_pp,j,i) - u(k-1,j,i) ) &
[888]	3075	+ ( ibit17 * adv_mom_5 &
	3076	) * &
[861]	3077	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
[888]	3078	)
[713]	3079	!
[861]	3080	!-- Calculate the divergence of the velocity field. A respective
	3081	!-- correction is needed to overcome numerical instabilities caused
	3082	!-- by a not sufficient reduction of divergences near topography.
	3083	div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
	3084	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
	3085	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) &
	3086	* ddzw(k) &
	3087	) * 0.5
[713]	3088
[861]	3089	tend(k,j,i) = tend(k,j,i) - ( &
	3090	( flux_r(k) + diss_r(k) &
	3091	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
	3092	+ ( flux_n(k) + diss_n(k) &
	3093	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy &
	3094	+ ( flux_t(k) + diss_t(k) &
	3095	- flux_d - diss_d &
	3096	) * ddzw(k) &
	3097	) + div * u(k,j,i)
	3098
	3099	swap_flux_x_local_u(k,j) = flux_r(k)
	3100	swap_diss_x_local_u(k,j) = diss_r(k)
	3101	swap_flux_y_local_u(k) = flux_n(k)
	3102	swap_diss_y_local_u(k) = diss_n(k)
	3103	flux_d = flux_t(k)
	3104	diss_d = diss_t(k)
[713]	3105	!
[861]	3106	!-- Statistical Evaluation of u'u'. The factor has to be applied
	3107	!-- for right evaluation when gallilei_trans = .T. .
	3108	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
	3109	+ ( flux_r(k) * &
	3110	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	3111	/ ( u_comp(k) - gu + 1.0E-20 ) &
	3112	+ diss_r(k) * &
	3113	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
	3114	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
[743]	3115	* weight_substep(intermediate_timestep_count)
[861]	3116	!
	3117	!-- Statistical Evaluation of w'u'.
	3118	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
	3119	+ ( flux_t(k) + diss_t(k) ) &
	3120	* weight_substep(intermediate_timestep_count)
	3121	ENDDO
[713]	3122	ENDDO
	3123	ENDDO
[861]	3124	sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
[713]	3125
	3126
	3127	END SUBROUTINE advec_u_ws
	3128
	3129
	3130	!------------------------------------------------------------------------------!
[1015]	3131	! Advection of u - Call for all grid points - accelerator version
	3132	!------------------------------------------------------------------------------!
	3133	SUBROUTINE advec_u_ws_acc
	3134
	3135	USE arrays_3d
	3136	USE constants
	3137	USE control_parameters
	3138	USE grid_variables
	3139	USE indices
	3140	USE statistics
	3141
	3142	IMPLICIT NONE
	3143
	3144	INTEGER :: i, ibit9, ibit10, ibit11, ibit12, ibit13, ibit14, ibit15, &
	3145	ibit16, ibit17, j, k, k_mmm, k_mm, k_pp, k_ppp, tn = 0
	3146
	3147	REAL :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, &
	3148	flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
	3149	u_comp, u_comp_l, v_comp, v_comp_s, w_comp
	3150
	3151
	3152	gu = 2.0 * u_gtrans
	3153	gv = 2.0 * v_gtrans
	3154
	3155	!
	3156	!-- Computation of fluxes and tendency terms
[1221]	3157	!$acc kernels present( ddzw, tend, u, v, w, wall_flags_0, wall_flags_00 )
[1015]	3158	!$acc loop
[1128]	3159	DO i = i_left, i_right
	3160	DO j = j_south, j_north
[1015]	3161	!$acc loop vector( 32 )
	3162	DO k = nzb+1, nzt
	3163
	3164	ibit11 = IBITS(wall_flags_0(k,j,i),11,1)
	3165	ibit10 = IBITS(wall_flags_0(k,j,i),10,1)
	3166	ibit9 = IBITS(wall_flags_0(k,j,i),9,1)
	3167
	3168	u_comp_l = u(k,j,i) + u(k,j,i-1) - gu
	3169	flux_l = u_comp_l * ( &
	3170	( 37.0 * ibit11 * adv_mom_5 &
	3171	+ 7.0 * ibit10 * adv_mom_3 &
	3172	+ ibit9 * adv_mom_1 &
	3173	) * &
	3174	( u(k,j,i) + u(k,j,i-1) ) &
	3175	- ( 8.0 * ibit11 * adv_mom_5 &
	3176	+ ibit10 * adv_mom_3 &
	3177	) * &
	3178	( u(k,j,i+1) + u(k,j,i-2) ) &
	3179	+ ( ibit11 * adv_mom_5 &
	3180	) * &
	3181	( u(k,j,i+2) + u(k,j,i-3) ) &
	3182	)
	3183
	3184	diss_l = - ABS( u_comp_l ) * ( &
	3185	( 10.0 * ibit11 * adv_mom_5 &
	3186	+ 3.0 * ibit10 * adv_mom_3 &
	3187	+ ibit9 * adv_mom_1 &
	3188	) * &
	3189	( u(k,j,i) - u(k,j,i-1) ) &
	3190	- ( 5.0 * ibit11 * adv_mom_5 &
	3191	+ ibit10 * adv_mom_3 &
	3192	) * &
	3193	( u(k,j,i+1) - u(k,j,i-2) ) &
	3194	+ ( ibit11 * adv_mom_5 &
	3195	) * &
	3196	( u(k,j,i+2) - u(k,j,i-3) ) &
	3197	)
	3198
	3199	u_comp = u(k,j,i+1) + u(k,j,i)
	3200	flux_r = ( u_comp - gu ) * ( &
	3201	( 37.0 * ibit11 * adv_mom_5 &
	3202	+ 7.0 * ibit10 * adv_mom_3 &
	3203	+ ibit9 * adv_mom_1 &
	3204	) * &
	3205	( u(k,j,i+1) + u(k,j,i) ) &
	3206	- ( 8.0 * ibit11 * adv_mom_5 &
	3207	+ ibit10 * adv_mom_3 &
	3208	) * &
	3209	( u(k,j,i+2) + u(k,j,i-1) ) &
	3210	+ ( ibit11 * adv_mom_5 &
	3211	) * &
	3212	( u(k,j,i+3) + u(k,j,i-2) ) &
	3213	)
	3214
	3215	diss_r = - ABS( u_comp - gu ) * ( &
	3216	( 10.0 * ibit11 * adv_mom_5 &
	3217	+ 3.0 * ibit10 * adv_mom_3 &
	3218	+ ibit9 * adv_mom_1 &
	3219	) * &
	3220	( u(k,j,i+1) - u(k,j,i) ) &
	3221	- ( 5.0 * ibit11 * adv_mom_5 &
	3222	+ ibit10 * adv_mom_3 &
	3223	) * &
	3224	( u(k,j,i+2) - u(k,j,i-1) ) &
	3225	+ ( ibit11 * adv_mom_5 &
	3226	) * &
	3227	( u(k,j,i+3) - u(k,j,i-2) ) &
	3228	)
	3229
	3230	ibit14 = IBITS(wall_flags_0(k,j,i),14,1)
	3231	ibit13 = IBITS(wall_flags_0(k,j,i),13,1)
	3232	ibit12 = IBITS(wall_flags_0(k,j,i),12,1)
	3233
	3234	v_comp_s = v(k,j,i) + v(k,j,i-1) - gv
	3235	flux_s = v_comp_s * ( &
	3236	( 37.0 * ibit14 * adv_mom_5 &
	3237	+ 7.0 * ibit13 * adv_mom_3 &
	3238	+ ibit12 * adv_mom_1 &
	3239	) * &
	3240	( u(k,j,i) + u(k,j-1,i) ) &
	3241	- ( 8.0 * ibit14 * adv_mom_5 &
	3242	+ ibit13 * adv_mom_3 &
	3243	) * &
	3244	( u(k,j+1,i) + u(k,j-2,i) ) &
	3245	+ ( ibit14 * adv_mom_5 &
	3246	) * &
	3247	( u(k,j+2,i) + u(k,j-3,i) ) &
	3248	)
	3249
	3250	diss_s = - ABS ( v_comp_s ) * ( &
	3251	( 10.0 * ibit14 * adv_mom_5 &
	3252	+ 3.0 * ibit13 * adv_mom_3 &
	3253	+ ibit12 * adv_mom_1 &
	3254	) * &
	3255	( u(k,j,i) - u(k,j-1,i) ) &
	3256	- ( 5.0 * ibit14 * adv_mom_5 &
	3257	+ ibit13 * adv_mom_3 &
	3258	) * &
	3259	( u(k,j+1,i) - u(k,j-2,i) ) &
	3260	+ ( ibit14 * adv_mom_5 &
	3261	) * &
	3262	( u(k,j+2,i) - u(k,j-3,i) ) &
	3263	)
	3264
	3265
	3266	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
	3267	flux_n = v_comp * ( &
	3268	( 37.0 * ibit14 * adv_mom_5 &
	3269	+ 7.0 * ibit13 * adv_mom_3 &
	3270	+ ibit12 * adv_mom_1 &
	3271	) * &
	3272	( u(k,j+1,i) + u(k,j,i) ) &
	3273	- ( 8.0 * ibit14 * adv_mom_5 &
	3274	+ ibit13 * adv_mom_3 &
	3275	) * &
	3276	( u(k,j+2,i) + u(k,j-1,i) ) &
	3277	+ ( ibit14 * adv_mom_5 &
	3278	) * &
	3279	( u(k,j+3,i) + u(k,j-2,i) ) &
	3280	)
	3281
	3282	diss_n = - ABS ( v_comp ) * ( &
	3283	( 10.0 * ibit14 * adv_mom_5 &
	3284	+ 3.0 * ibit13 * adv_mom_3 &
	3285	+ ibit12 * adv_mom_1 &
	3286	) * &
	3287	( u(k,j+1,i) - u(k,j,i) ) &
	3288	- ( 5.0 * ibit14 * adv_mom_5 &
	3289	+ ibit13 * adv_mom_3 &
	3290	) * &
	3291	( u(k,j+2,i) - u(k,j-1,i) ) &
	3292	+ ( ibit14 * adv_mom_5 &
	3293	) * &
	3294	( u(k,j+3,i) - u(k,j-2,i) ) &
	3295	)
	3296
	3297	ibit17 = IBITS(wall_flags_0(k-1,j,i),17,1)
	3298	ibit16 = IBITS(wall_flags_0(k-1,j,i),16,1)
	3299	ibit15 = IBITS(wall_flags_0(k-1,j,i),15,1)
	3300
[1027]	3301	k_pp = k + 2 * ibit17
	3302	k_mm = k - 2 * ( ibit16 + ibit17 )
[1015]	3303	k_mmm = k - 3 * ibit17
	3304
	3305	w_comp = w(k-1,j,i) + w(k-1,j,i-1)
	3306	flux_d = w_comp * ( &
	3307	( 37.0 * ibit17 * adv_mom_5 &
	3308	+ 7.0 * ibit16 * adv_mom_3 &
	3309	+ ibit15 * adv_mom_1 &
	3310	) * &
	3311	( u(k,j,i) + u(k-1,j,i) ) &
	3312	- ( 8.0 * ibit17 * adv_mom_5 &
	3313	+ ibit16 * adv_mom_3 &
	3314	) * &
	3315	( u(k+1,j,i) + u(k_mm,j,i) ) &
	3316	+ ( ibit17 * adv_mom_5 &
	3317	) * &
	3318	( u(k_pp,j,i) + u(k_mmm,j,i) ) &
	3319	)
	3320
	3321	diss_d = - ABS( w_comp ) * ( &
	3322	( 10.0 * ibit17 * adv_mom_5 &
	3323	+ 3.0 * ibit16 * adv_mom_3 &
	3324	+ ibit15 * adv_mom_1 &
	3325	) * &
	3326	( u(k,j,i) - u(k-1,j,i) ) &
	3327	- ( 5.0 * ibit17 * adv_mom_5 &
	3328	+ ibit16 * adv_mom_3 &
	3329	) * &
	3330	( u(k+1,j,i) - u(k_mm,j,i) ) &
	3331	+ ( ibit17 * adv_mom_5 &
	3332	) * &
	3333	( u(k_pp,j,i) - u(k_mmm,j,i) ) &
	3334	)
	3335	!
	3336	!-- k index has to be modified near bottom and top, else array
	3337	!-- subscripts will be exceeded.
	3338	ibit17 = IBITS(wall_flags_0(k,j,i),17,1)
	3339	ibit16 = IBITS(wall_flags_0(k,j,i),16,1)
	3340	ibit15 = IBITS(wall_flags_0(k,j,i),15,1)
	3341
	3342	k_ppp = k + 3 * ibit17
	3343	k_pp = k + 2 * ( 1 - ibit15 )
	3344	k_mm = k - 2 * ibit17
	3345
	3346	w_comp = w(k,j,i) + w(k,j,i-1)
	3347	flux_t = w_comp * ( &
	3348	( 37.0 * ibit17 * adv_mom_5 &
	3349	+ 7.0 * ibit16 * adv_mom_3 &
	3350	+ ibit15 * adv_mom_1 &
	3351	) * &
	3352	( u(k+1,j,i) + u(k,j,i) ) &
	3353	- ( 8.0 * ibit17 * adv_mom_5 &
	3354	+ ibit16 * adv_mom_3 &
	3355	) * &
	3356	( u(k_pp,j,i) + u(k-1,j,i) ) &
	3357	+ ( ibit17 * adv_mom_5 &
	3358	) * &
	3359	( u(k_ppp,j,i) + u(k_mm,j,i) ) &
	3360	)
	3361
	3362	diss_t = - ABS( w_comp ) * ( &
	3363	( 10.0 * ibit17 * adv_mom_5 &
	3364	+ 3.0 * ibit16 * adv_mom_3 &
	3365	+ ibit15 * adv_mom_1 &
	3366	) * &
	3367	( u(k+1,j,i) - u(k,j,i) ) &
	3368	- ( 5.0 * ibit17 * adv_mom_5 &
	3369	+ ibit16 * adv_mom_3 &
	3370	) * &
	3371	( u(k_pp,j,i) - u(k-1,j,i) ) &
	3372	+ ( ibit17 * adv_mom_5 &
	3373	) * &
	3374	( u(k_ppp,j,i) - u(k_mm,j,i) ) &
	3375	)
	3376	!
	3377	!-- Calculate the divergence of the velocity field. A respective
	3378	!-- correction is needed to overcome numerical instabilities caused
	3379	!-- by a not sufficient reduction of divergences near topography.
	3380	div = ( ( u_comp - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &
	3381	+ ( v_comp + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &
	3382	+ ( w_comp - ( w(k-1,j,i) + w(k-1,j,i-1) ) ) &
	3383	* ddzw(k) &
	3384	) * 0.5
	3385
	3386	tend(k,j,i) = - ( &
	3387	( flux_r + diss_r - flux_l - diss_l ) * ddx &
	3388	+ ( flux_n + diss_n - flux_s - diss_s ) * ddy &
	3389	+ ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
	3390	) + div * u(k,j,i)
	3391
	3392	!++
	3393	!-- Statistical Evaluation of u'u'. The factor has to be applied
	3394	!-- for right evaluation when gallilei_trans = .T. .
	3395	! sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
	3396	! + ( flux_r * &
	3397	! ( u_comp - 2.0 * hom(k,1,1,0) ) &
	3398	! / ( u_comp - gu + 1.0E-20 ) &
	3399	! + diss_r * &
	3400	! ABS( u_comp - 2.0 * hom(k,1,1,0) ) &
	3401	! / ( ABS( u_comp - gu ) + 1.0E-20 ) ) &
	3402	! * weight_substep(intermediate_timestep_count)
	3403	!
	3404	!-- Statistical Evaluation of w'u'.
	3405	! sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
	3406	! + ( flux_t + diss_t ) &
	3407	! * weight_substep(intermediate_timestep_count)
	3408	ENDDO
	3409	ENDDO
	3410	ENDDO
	3411	!$acc end kernels
	3412
	3413	!++
	3414	! sums_us2_ws_l(nzb,tn) = sums_us2_ws_l(nzb+1,tn)
	3415
	3416	END SUBROUTINE advec_u_ws_acc
	3417
	3418
	3419	!------------------------------------------------------------------------------!
[713]	3420	! Advection of v - Call for all grid points
	3421	!------------------------------------------------------------------------------!
	3422	SUBROUTINE advec_v_ws
	3423
	3424	USE arrays_3d
	3425	USE constants
	3426	USE control_parameters
	3427	USE grid_variables
	3428	USE indices
	3429	USE statistics
	3430
	3431	IMPLICIT NONE
	3432
	3433
[888]	3434	INTEGER :: i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
	3435	ibit25, ibit26, j, k, k_mm, k_pp, k_ppp, tn = 0
[862]	3436	REAL :: diss_d, div, flux_d, gu, gv, u_comp, w_comp
	3437	REAL, DIMENSION(nzb+1:nzt) :: swap_diss_y_local_v, swap_flux_y_local_v
	3438	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_v, &
	3439	swap_flux_x_local_v
	3440	REAL, DIMENSION(nzb:nzt) :: diss_n, diss_r, diss_t, flux_n, flux_r, &
	3441	flux_t, v_comp
[713]	3442
	3443	gu = 2.0 * u_gtrans
	3444	gv = 2.0 * v_gtrans
	3445	!
[861]	3446	!-- First compute the whole left boundary of the processor domain
[713]	3447	i = nxl
	3448	DO j = nysv, nyn
[861]	3449	DO k = nzb+1, nzb_max
	3450
[888]	3451	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
	3452	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
	3453	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
	3454
[861]	3455	u_comp = u(k,j-1,i) + u(k,j,i) - gu
	3456	swap_flux_x_local_v(k,j) = u_comp * ( &
[888]	3457	( 37.0 * ibit20 * adv_mom_5 &
	3458	+ 7.0 * ibit19 * adv_mom_3 &
	3459	+ ibit18 * adv_mom_1 &
	3460	) * &
[861]	3461	( v(k,j,i) + v(k,j,i-1) ) &
[888]	3462	- ( 8.0 * ibit20 * adv_mom_5 &
	3463	+ ibit19 * adv_mom_3 &
	3464	) * &
[861]	3465	( v(k,j,i+1) + v(k,j,i-2) ) &
[888]	3466	+ ( ibit20 * adv_mom_5 &
	3467	) * &
[861]	3468	( v(k,j,i+2) + v(k,j,i-3) ) &
	3469	)
	3470
	3471	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
[888]	3472	( 10.0 * ibit20 * adv_mom_5 &
	3473	+ 3.0 * ibit19 * adv_mom_3 &
	3474	+ ibit18 * adv_mom_1 &
	3475	) * &
[861]	3476	( v(k,j,i) - v(k,j,i-1) ) &
[888]	3477	- ( 5.0 * ibit20 * adv_mom_5 &
	3478	+ ibit19 * adv_mom_3 &
	3479	) * &
[861]	3480	( v(k,j,i+1) - v(k,j,i-2) ) &
[888]	3481	+ ( ibit20 * adv_mom_5 &
	3482	) * &
[861]	3483	( v(k,j,i+2) - v(k,j,i-3) ) &
	3484	)
	3485
	3486	ENDDO
	3487
	3488	DO k = nzb_max+1, nzt
	3489
	3490	u_comp = u(k,j-1,i) + u(k,j,i) - gu
	3491	swap_flux_x_local_v(k,j) = u_comp * ( &
	3492	37.0 * ( v(k,j,i) + v(k,j,i-1) ) &
	3493	- 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) &
	3494	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
	3495	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
	3496	10.0 * ( v(k,j,i) - v(k,j,i-1) ) &
	3497	- 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) &
	3498	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
	3499
	3500	ENDDO
	3501
[713]	3502	ENDDO
	3503
	3504	DO i = nxl, nxr
	3505
[861]	3506	j = nysv
	3507	DO k = nzb+1, nzb_max
[713]	3508
[888]	3509	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
	3510	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
	3511	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
	3512
[861]	3513	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
	3514	swap_flux_y_local_v(k) = v_comp(k) * ( &
[888]	3515	( 37.0 * ibit23 * adv_mom_5 &
	3516	+ 7.0 * ibit22 * adv_mom_3 &
	3517	+ ibit21 * adv_mom_1 &
	3518	) * &
[861]	3519	( v(k,j,i) + v(k,j-1,i) ) &
[888]	3520	- ( 8.0 * ibit23 * adv_mom_5 &
	3521	+ ibit22 * adv_mom_3 &
	3522	) * &
[861]	3523	( v(k,j+1,i) + v(k,j-2,i) ) &
[888]	3524	+ ( ibit23 * adv_mom_5 &
	3525	) * &
[861]	3526	( v(k,j+2,i) + v(k,j-3,i) ) &
	3527	)
[713]	3528
[861]	3529	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
[888]	3530	( 10.0 * ibit23 * adv_mom_5 &
	3531	+ 3.0 * ibit22 * adv_mom_3 &
	3532	+ ibit21 * adv_mom_1 &
	3533	) * &
[861]	3534	( v(k,j,i) - v(k,j-1,i) ) &
[888]	3535	- ( 5.0 * ibit23 * adv_mom_5 &
	3536	+ ibit22 * adv_mom_3 &
	3537	) * &
[861]	3538	( v(k,j+1,i) - v(k,j-2,i) ) &
[888]	3539	+ ( ibit23 * adv_mom_5 &
	3540	) * &
[861]	3541	( v(k,j+2,i) - v(k,j-3,i) ) &
	3542	)
	3543
	3544	ENDDO
	3545
	3546	DO k = nzb_max+1, nzt
	3547
	3548	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
	3549	swap_flux_y_local_v(k) = v_comp(k) * ( &
	3550	37.0 * ( v(k,j,i) + v(k,j-1,i) ) &
	3551	- 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) &
	3552	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
	3553	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
	3554	10.0 * ( v(k,j,i) - v(k,j-1,i) ) &
	3555	- 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) &
	3556	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
	3557
	3558	ENDDO
	3559
[713]	3560	DO j = nysv, nyn
[861]	3561
	3562	flux_t(0) = 0.0
	3563	diss_t(0) = 0.0
	3564	flux_d = 0.0
	3565	diss_d = 0.0
	3566
	3567	DO k = nzb+1, nzb_max
	3568
[888]	3569	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
	3570	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
	3571	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
	3572
[861]	3573	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
	3574	flux_r(k) = u_comp * ( &
[888]	3575	( 37.0 * ibit20 * adv_mom_5 &
	3576	+ 7.0 * ibit19 * adv_mom_3 &
	3577	+ ibit18 * adv_mom_1 &
	3578	) * &
[861]	3579	( v(k,j,i+1) + v(k,j,i) ) &
[888]	3580	- ( 8.0 * ibit20 * adv_mom_5 &
	3581	+ ibit19 * adv_mom_3 &
	3582	) * &
[861]	3583	( v(k,j,i+2) + v(k,j,i-1) ) &
[888]	3584	+ ( ibit20 * adv_mom_5 &
	3585	) * &
[861]	3586	( v(k,j,i+3) + v(k,j,i-2) ) &
	3587	)
	3588
	3589	diss_r(k) = - ABS( u_comp ) * ( &
[888]	3590	( 10.0 * ibit20 * adv_mom_5 &
	3591	+ 3.0 * ibit19 * adv_mom_3 &
	3592	+ ibit18 * adv_mom_1 &
	3593	) * &
[861]	3594	( v(k,j,i+1) - v(k,j,i) ) &
[888]	3595	- ( 5.0 * ibit20 * adv_mom_5 &
	3596	+ ibit19 * adv_mom_3 &
	3597	) * &
[861]	3598	( v(k,j,i+2) - v(k,j,i-1) ) &
[888]	3599	+ ( ibit20 * adv_mom_5 &
	3600	) * &
[861]	3601	( v(k,j,i+3) - v(k,j,i-2) ) &
	3602	)
	3603
[888]	3604	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
	3605	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
	3606	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
	3607
[861]	3608	v_comp(k) = v(k,j+1,i) + v(k,j,i)
	3609	flux_n(k) = ( v_comp(k) - gv ) * ( &
[888]	3610	( 37.0 * ibit23 * adv_mom_5 &
	3611	+ 7.0 * ibit22 * adv_mom_3 &
	3612	+ ibit21 * adv_mom_1 &
	3613	) * &
[861]	3614	( v(k,j+1,i) + v(k,j,i) ) &
[888]	3615	- ( 8.0 * ibit23 * adv_mom_5 &
	3616	+ ibit22 * adv_mom_3 &
	3617	) * &
[861]	3618	( v(k,j+2,i) + v(k,j-1,i) ) &
[888]	3619	+ ( ibit23 * adv_mom_5 &
	3620	) * &
[861]	3621	( v(k,j+3,i) + v(k,j-2,i) ) &
	3622	)
	3623
	3624	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
[888]	3625	( 10.0 * ibit23 * adv_mom_5 &
	3626	+ 3.0 * ibit22 * adv_mom_3 &
	3627	+ ibit21 * adv_mom_1 &
	3628	) * &
[861]	3629	( v(k,j+1,i) - v(k,j,i) ) &
[888]	3630	- ( 5.0 * ibit23 * adv_mom_5 &
	3631	+ ibit22 * adv_mom_3 &
	3632	) * &
[861]	3633	( v(k,j+2,i) - v(k,j-1,i) ) &
[888]	3634	+ ( ibit23 * adv_mom_5 &
	3635	) * &
[861]	3636	( v(k,j+3,i) - v(k,j-2,i) ) &
	3637	)
[713]	3638	!
[861]	3639	!-- k index has to be modified near bottom and top, else array
	3640	!-- subscripts will be exceeded.
[888]	3641	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
	3642	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
	3643	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
[713]	3644
[888]	3645	k_ppp = k + 3 * ibit26
	3646	k_pp = k + 2 * ( 1 - ibit24 )
	3647	k_mm = k - 2 * ibit26
	3648
[861]	3649	w_comp = w(k,j-1,i) + w(k,j,i)
	3650	flux_t(k) = w_comp * ( &
[888]	3651	( 37.0 * ibit26 * adv_mom_5 &
	3652	+ 7.0 * ibit25 * adv_mom_3 &
	3653	+ ibit24 * adv_mom_1 &
	3654	) * &
[861]	3655	( v(k+1,j,i) + v(k,j,i) ) &
[888]	3656	- ( 8.0 * ibit26 * adv_mom_5 &
	3657	+ ibit25 * adv_mom_3 &
	3658	) * &
[861]	3659	( v(k_pp,j,i) + v(k-1,j,i) ) &
[888]	3660	+ ( ibit26 * adv_mom_5 &
	3661	) * &
[861]	3662	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
	3663	)
	3664
	3665	diss_t(k) = - ABS( w_comp ) * ( &
[888]	3666	( 10.0 * ibit26 * adv_mom_5 &
	3667	+ 3.0 * ibit25 * adv_mom_3 &
	3668	+ ibit24 * adv_mom_1 &
	3669	) * &
[861]	3670	( v(k+1,j,i) - v(k,j,i) ) &
[888]	3671	- ( 5.0 * ibit26 * adv_mom_5 &
	3672	+ ibit25 * adv_mom_3 &
	3673	) * &
[861]	3674	( v(k_pp,j,i) - v(k-1,j,i) ) &
[888]	3675	+ ( ibit26 * adv_mom_5 &
	3676	) * &
[861]	3677	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
	3678	)
[713]	3679	!
[861]	3680	!-- Calculate the divergence of the velocity field. A respective
	3681	!-- correction is needed to overcome numerical instabilities caused
	3682	!-- by a not sufficient reduction of divergences near topography.
	3683	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
	3684	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
	3685	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) &
	3686	) * ddzw(k) &
	3687	) * 0.5
[713]	3688
[861]	3689	tend(k,j,i) = tend(k,j,i) - ( &
	3690	( flux_r(k) + diss_r(k) &
	3691	- swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) &
	3692	) * ddx &
	3693	+ ( flux_n(k) + diss_n(k) &
	3694	- swap_flux_y_local_v(k) - swap_diss_y_local_v(k) &
	3695	) * ddy &
	3696	+ ( flux_t(k) + diss_t(k) &
	3697	- flux_d - diss_d &
	3698	) * ddzw(k) &
	3699	) + v(k,j,i) * div
	3700
	3701	swap_flux_x_local_v(k,j) = flux_r(k)
	3702	swap_diss_x_local_v(k,j) = diss_r(k)
	3703	swap_flux_y_local_v(k) = flux_n(k)
	3704	swap_diss_y_local_v(k) = diss_n(k)
	3705	flux_d = flux_t(k)
	3706	diss_d = diss_t(k)
	3707
[713]	3708	!
[861]	3709	!-- Statistical Evaluation of v'v'. The factor has to be applied
	3710	!-- for right evaluation when gallilei_trans = .T. .
	3711	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
	3712	+ ( flux_n(k) &
	3713	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
	3714	/ ( v_comp(k) - gv + 1.0E-20 ) &
	3715	+ diss_n(k) &
	3716	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
	3717	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
	3718	* weight_substep(intermediate_timestep_count)
[713]	3719	!
[861]	3720	!-- Statistical Evaluation of w'v'.
	3721	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
	3722	+ ( flux_t(k) + diss_t(k) ) &
	3723	* weight_substep(intermediate_timestep_count)
	3724
	3725	ENDDO
	3726
	3727	DO k = nzb_max+1, nzt
	3728
	3729	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
	3730	flux_r(k) = u_comp * ( &
	3731	37.0 * ( v(k,j,i+1) + v(k,j,i) ) &
	3732	- 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) &
	3733	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
	3734
	3735	diss_r(k) = - ABS( u_comp ) * ( &
	3736	10.0 * ( v(k,j,i+1) - v(k,j,i) ) &
	3737	- 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) &
	3738	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
	3739
	3740
	3741	v_comp(k) = v(k,j+1,i) + v(k,j,i)
	3742	flux_n(k) = ( v_comp(k) - gv ) * ( &
	3743	37.0 * ( v(k,j+1,i) + v(k,j,i) ) &
	3744	- 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) &
	3745	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
	3746
	3747	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
	3748	10.0 * ( v(k,j+1,i) - v(k,j,i) ) &
	3749	- 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) &
	3750	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
[713]	3751	!
[861]	3752	!-- k index has to be modified near bottom and top, else array
	3753	!-- subscripts will be exceeded.
[888]	3754	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
	3755	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
	3756	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
[713]	3757
[888]	3758	k_ppp = k + 3 * ibit26
	3759	k_pp = k + 2 * ( 1 - ibit24 )
	3760	k_mm = k - 2 * ibit26
	3761
[861]	3762	w_comp = w(k,j-1,i) + w(k,j,i)
	3763	flux_t(k) = w_comp * ( &
[888]	3764	( 37.0 * ibit26 * adv_mom_5 &
	3765	+ 7.0 * ibit25 * adv_mom_3 &
	3766	+ ibit24 * adv_mom_1 &
	3767	) * &
[861]	3768	( v(k+1,j,i) + v(k,j,i) ) &
[888]	3769	- ( 8.0 * ibit26 * adv_mom_5 &
	3770	+ ibit25 * adv_mom_3 &
	3771	) * &
[861]	3772	( v(k_pp,j,i) + v(k-1,j,i) ) &
[888]	3773	+ ( ibit26 * adv_mom_5 &
	3774	) * &
[861]	3775	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
[888]	3776	)
[861]	3777
	3778	diss_t(k) = - ABS( w_comp ) * ( &
[888]	3779	( 10.0 * ibit26 * adv_mom_5 &
	3780	+ 3.0 * ibit25 * adv_mom_3 &
	3781	+ ibit24 * adv_mom_1 &
	3782	) * &
[861]	3783	( v(k+1,j,i) - v(k,j,i) ) &
[888]	3784	- ( 5.0 * ibit26 * adv_mom_5 &
	3785	+ ibit25 * adv_mom_3 &
	3786	) * &
[861]	3787	( v(k_pp,j,i) - v(k-1,j,i) ) &
[888]	3788	+ ( ibit26 * adv_mom_5 &
	3789	) * &
[861]	3790	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
[888]	3791	)
[713]	3792	!
[861]	3793	!-- Calculate the divergence of the velocity field. A respective
	3794	!-- correction is needed to overcome numerical instabilities caused
	3795	!-- by a not sufficient reduction of divergences near topography.
	3796	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
	3797	+ ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
	3798	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) ) &
	3799	* ddzw(k) &
	3800	) * 0.5
	3801
	3802	tend(k,j,i) = tend(k,j,i) - ( &
	3803	( flux_r(k) + diss_r(k) &
	3804	- swap_flux_x_local_v(k,j) - swap_diss_x_local_v(k,j) &
	3805	) * ddx &
	3806	+ ( flux_n(k) + diss_n(k) &
	3807	- swap_flux_y_local_v(k) - swap_diss_y_local_v(k) &
	3808	) * ddy &
	3809	+ ( flux_t(k) + diss_t(k) &
	3810	- flux_d - diss_d &
	3811	) * ddzw(k) &
	3812	) + v(k,j,i) * div
	3813
	3814	swap_flux_x_local_v(k,j) = flux_r(k)
	3815	swap_diss_x_local_v(k,j) = diss_r(k)
	3816	swap_flux_y_local_v(k) = flux_n(k)
	3817	swap_diss_y_local_v(k) = diss_n(k)
	3818	flux_d = flux_t(k)
	3819	diss_d = diss_t(k)
	3820
	3821	!
	3822	!-- Statistical Evaluation of v'v'. The factor has to be applied
	3823	!-- for right evaluation when gallilei_trans = .T. .
	3824	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
	3825	+ ( flux_n(k) &
	3826	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
	3827	/ ( v_comp(k) - gv + 1.0E-20 ) &
	3828	+ diss_n(k) &
	3829	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
	3830	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
	3831	* weight_substep(intermediate_timestep_count)
	3832	!
	3833	!-- Statistical Evaluation of w'v'.
	3834	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
	3835	+ ( flux_t(k) + diss_t(k) ) &
	3836	* weight_substep(intermediate_timestep_count)
	3837
	3838	ENDDO
[713]	3839	ENDDO
	3840	ENDDO
[861]	3841	sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
[713]	3842
[861]	3843
[713]	3844	END SUBROUTINE advec_v_ws
	3845
	3846
	3847	!------------------------------------------------------------------------------!
[1015]	3848	! Advection of v - Call for all grid points - accelerator version
	3849	!------------------------------------------------------------------------------!
	3850	SUBROUTINE advec_v_ws_acc
	3851
	3852	USE arrays_3d
	3853	USE constants
	3854	USE control_parameters
	3855	USE grid_variables
	3856	USE indices
	3857	USE statistics
	3858
	3859	IMPLICIT NONE
	3860
	3861
	3862	INTEGER :: i, ibit18, ibit19, ibit20, ibit21, ibit22, ibit23, ibit24, &
	3863	ibit25, ibit26, j, k, k_mm, k_mmm, k_pp, k_ppp, tn = 0
	3864
	3865	REAL :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, &
	3866	flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
	3867	u_comp, u_comp_l, v_comp, v_comp_s, w_comp
	3868
	3869	gu = 2.0 * u_gtrans
	3870	gv = 2.0 * v_gtrans
	3871
	3872	!
	3873	!-- Computation of fluxes and tendency terms
[1221]	3874	!$acc kernels present( ddzw, tend, u, v, w, wall_flags_0, wall_flags_00 )
[1015]	3875	!$acc loop
[1128]	3876	DO i = i_left, i_right
	3877	DO j = j_south, j_north
[1015]	3878	!$acc loop vector( 32 )
	3879	DO k = nzb+1, nzt
	3880
	3881	ibit20 = IBITS(wall_flags_0(k,j,i),20,1)
	3882	ibit19 = IBITS(wall_flags_0(k,j,i),19,1)
	3883	ibit18 = IBITS(wall_flags_0(k,j,i),18,1)
	3884
	3885	u_comp_l = u(k,j-1,i) + u(k,j,i) - gu
	3886	flux_l = u_comp_l * ( &
	3887	( 37.0 * ibit20 * adv_mom_5 &
	3888	+ 7.0 * ibit19 * adv_mom_3 &
	3889	+ ibit18 * adv_mom_1 &
	3890	) * &
	3891	( v(k,j,i) + v(k,j,i-1) ) &
	3892	- ( 8.0 * ibit20 * adv_mom_5 &
	3893	+ ibit19 * adv_mom_3 &
	3894	) * &
	3895	( v(k,j,i+1) + v(k,j,i-2) ) &
	3896	+ ( ibit20 * adv_mom_5 &
	3897	) * &
	3898	( v(k,j,i+2) + v(k,j,i-3) ) &
	3899	)
	3900
	3901	diss_l = - ABS( u_comp_l ) * ( &
	3902	( 10.0 * ibit20 * adv_mom_5 &
	3903	+ 3.0 * ibit19 * adv_mom_3 &
	3904	+ ibit18 * adv_mom_1 &
	3905	) * &
	3906	( v(k,j,i) - v(k,j,i-1) ) &
	3907	- ( 5.0 * ibit20 * adv_mom_5 &
	3908	+ ibit19 * adv_mom_3 &
	3909	) * &
	3910	( v(k,j,i+1) - v(k,j,i-2) ) &
	3911	+ ( ibit20 * adv_mom_5 &
	3912	) * &
	3913	( v(k,j,i+2) - v(k,j,i-3) ) &
	3914	)
	3915
	3916	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
	3917	flux_r = u_comp * ( &
	3918	( 37.0 * ibit20 * adv_mom_5 &
	3919	+ 7.0 * ibit19 * adv_mom_3 &
	3920	+ ibit18 * adv_mom_1 &
	3921	) * &
	3922	( v(k,j,i+1) + v(k,j,i) ) &
	3923	- ( 8.0 * ibit20 * adv_mom_5 &
	3924	+ ibit19 * adv_mom_3 &
	3925	) * &
	3926	( v(k,j,i+2) + v(k,j,i-1) ) &
	3927	+ ( ibit20 * adv_mom_5 &
	3928	) * &
	3929	( v(k,j,i+3) + v(k,j,i-2) ) &
	3930	)
	3931
	3932	diss_r = - ABS( u_comp ) * ( &
	3933	( 10.0 * ibit20 * adv_mom_5 &
	3934	+ 3.0 * ibit19 * adv_mom_3 &
	3935	+ ibit18 * adv_mom_1 &
	3936	) * &
	3937	( v(k,j,i+1) - v(k,j,i) ) &
	3938	- ( 5.0 * ibit20 * adv_mom_5 &
	3939	+ ibit19 * adv_mom_3 &
	3940	) * &
	3941	( v(k,j,i+2) - v(k,j,i-1) ) &
	3942	+ ( ibit20 * adv_mom_5 &
	3943	) * &
	3944	( v(k,j,i+3) - v(k,j,i-2) ) &
	3945	)
	3946
	3947	ibit23 = IBITS(wall_flags_0(k,j,i),23,1)
	3948	ibit22 = IBITS(wall_flags_0(k,j,i),22,1)
	3949	ibit21 = IBITS(wall_flags_0(k,j,i),21,1)
	3950
	3951
	3952	v_comp_s = v(k,j,i) + v(k,j-1,i) - gv
	3953	flux_s = v_comp_s * ( &
	3954	( 37.0 * ibit23 * adv_mom_5 &
	3955	+ 7.0 * ibit22 * adv_mom_3 &
	3956	+ ibit21 * adv_mom_1 &
	3957	) * &
	3958	( v(k,j,i) + v(k,j-1,i) ) &
	3959	- ( 8.0 * ibit23 * adv_mom_5 &
	3960	+ ibit22 * adv_mom_3 &
	3961	) * &
	3962	( v(k,j+1,i) + v(k,j-2,i) ) &
	3963	+ ( ibit23 * adv_mom_5 &
	3964	) * &
	3965	( v(k,j+2,i) + v(k,j-3,i) ) &
	3966	)
	3967
	3968	diss_s = - ABS( v_comp_s ) * ( &
	3969	( 10.0 * ibit23 * adv_mom_5 &
	3970	+ 3.0 * ibit22 * adv_mom_3 &
	3971	+ ibit21 * adv_mom_1 &
	3972	) * &
	3973	( v(k,j,i) - v(k,j-1,i) ) &
	3974	- ( 5.0 * ibit23 * adv_mom_5 &
	3975	+ ibit22 * adv_mom_3 &
	3976	) * &
	3977	( v(k,j+1,i) - v(k,j-2,i) ) &
	3978	+ ( ibit23 * adv_mom_5 &
	3979	) * &
	3980	( v(k,j+2,i) - v(k,j-3,i) ) &
	3981	)
	3982
	3983	v_comp = v(k,j+1,i) + v(k,j,i)
	3984	flux_n = ( v_comp - gv ) * ( &
	3985	( 37.0 * ibit23 * adv_mom_5 &
	3986	+ 7.0 * ibit22 * adv_mom_3 &
	3987	+ ibit21 * adv_mom_1 &
	3988	) * &
	3989	( v(k,j+1,i) + v(k,j,i) ) &
	3990	- ( 8.0 * ibit23 * adv_mom_5 &
	3991	+ ibit22 * adv_mom_3 &
	3992	) * &
	3993	( v(k,j+2,i) + v(k,j-1,i) ) &
	3994	+ ( ibit23 * adv_mom_5 &
	3995	) * &
	3996	( v(k,j+3,i) + v(k,j-2,i) ) &
	3997	)
	3998
	3999	diss_n = - ABS( v_comp - gv ) * ( &
	4000	( 10.0 * ibit23 * adv_mom_5 &
	4001	+ 3.0 * ibit22 * adv_mom_3 &
	4002	+ ibit21 * adv_mom_1 &
	4003	) * &
	4004	( v(k,j+1,i) - v(k,j,i) ) &
	4005	- ( 5.0 * ibit23 * adv_mom_5 &
	4006	+ ibit22 * adv_mom_3 &
	4007	) * &
	4008	( v(k,j+2,i) - v(k,j-1,i) ) &
	4009	+ ( ibit23 * adv_mom_5 &
	4010	) * &
	4011	( v(k,j+3,i) - v(k,j-2,i) ) &
	4012	)
	4013
	4014	ibit26 = IBITS(wall_flags_0(k-1,j,i),26,1)
	4015	ibit25 = IBITS(wall_flags_0(k-1,j,i),25,1)
	4016	ibit24 = IBITS(wall_flags_0(k-1,j,i),24,1)
	4017
[1027]	4018	k_pp = k + 2 * ibit26
	4019	k_mm = k - 2 * ( ibit25 + ibit26 )
[1015]	4020	k_mmm = k - 3 * ibit26
	4021
	4022	w_comp = w(k-1,j-1,i) + w(k-1,j,i)
	4023	flux_d = w_comp * ( &
	4024	( 37.0 * ibit26 * adv_mom_5 &
	4025	+ 7.0 * ibit25 * adv_mom_3 &
	4026	+ ibit24 * adv_mom_1 &
	4027	) * &
	4028	( v(k,j,i) + v(k-1,j,i) ) &
	4029	- ( 8.0 * ibit26 * adv_mom_5 &
	4030	+ ibit25 * adv_mom_3 &
	4031	) * &
	4032	( v(k+1,j,i) + v(k_mm,j,i) ) &
	4033	+ ( ibit26 * adv_mom_5 &
	4034	) * &
	4035	( v(k_pp,j,i) + v(k_mmm,j,i) ) &
	4036	)
	4037
	4038	diss_d = - ABS( w_comp ) * ( &
	4039	( 10.0 * ibit26 * adv_mom_5 &
	4040	+ 3.0 * ibit25 * adv_mom_3 &
	4041	+ ibit24 * adv_mom_1 &
	4042	) * &
	4043	( v(k,j,i) - v(k-1,j,i) ) &
	4044	- ( 5.0 * ibit26 * adv_mom_5 &
	4045	+ ibit25 * adv_mom_3 &
	4046	) * &
	4047	( v(k+1,j,i) - v(k_mm,j,i) ) &
	4048	+ ( ibit26 * adv_mom_5 &
	4049	) * &
	4050	( v(k_pp,j,i) - v(k_mmm,j,i) ) &
	4051	)
	4052	!
	4053	!-- k index has to be modified near bottom and top, else array
	4054	!-- subscripts will be exceeded.
	4055	ibit26 = IBITS(wall_flags_0(k,j,i),26,1)
	4056	ibit25 = IBITS(wall_flags_0(k,j,i),25,1)
	4057	ibit24 = IBITS(wall_flags_0(k,j,i),24,1)
	4058
	4059	k_ppp = k + 3 * ibit26
	4060	k_pp = k + 2 * ( 1 - ibit24 )
	4061	k_mm = k - 2 * ibit26
	4062
	4063	w_comp = w(k,j-1,i) + w(k,j,i)
	4064	flux_t = w_comp * ( &
	4065	( 37.0 * ibit26 * adv_mom_5 &
	4066	+ 7.0 * ibit25 * adv_mom_3 &
	4067	+ ibit24 * adv_mom_1 &
	4068	) * &
	4069	( v(k+1,j,i) + v(k,j,i) ) &
	4070	- ( 8.0 * ibit26 * adv_mom_5 &
	4071	+ ibit25 * adv_mom_3 &
	4072	) * &
	4073	( v(k_pp,j,i) + v(k-1,j,i) ) &
	4074	+ ( ibit26 * adv_mom_5 &
	4075	) * &
	4076	( v(k_ppp,j,i) + v(k_mm,j,i) ) &
	4077	)
	4078
	4079	diss_t = - ABS( w_comp ) * ( &
	4080	( 10.0 * ibit26 * adv_mom_5 &
	4081	+ 3.0 * ibit25 * adv_mom_3 &
	4082	+ ibit24 * adv_mom_1 &
	4083	) * &
	4084	( v(k+1,j,i) - v(k,j,i) ) &
	4085	- ( 5.0 * ibit26 * adv_mom_5 &
	4086	+ ibit25 * adv_mom_3 &
	4087	) * &
	4088	( v(k_pp,j,i) - v(k-1,j,i) ) &
	4089	+ ( ibit26 * adv_mom_5 &
	4090	) * &
	4091	( v(k_ppp,j,i) - v(k_mm,j,i) ) &
	4092	)
	4093	!
	4094	!-- Calculate the divergence of the velocity field. A respective
	4095	!-- correction is needed to overcome numerical instabilities caused
	4096	!-- by a not sufficient reduction of divergences near topography.
	4097	div = ( ( u_comp + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &
	4098	+ ( v_comp - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &
	4099	+ ( w_comp - ( w(k-1,j-1,i) + w(k-1,j,i) ) &
	4100	) * ddzw(k) &
	4101	) * 0.5
	4102
	4103	tend(k,j,i) = - ( &
	4104	( flux_r + diss_r - flux_l - diss_l ) * ddx &
	4105	+ ( flux_n + diss_n - flux_s - diss_s ) * ddy &
	4106	+ ( flux_t + diss_t - flux_d - diss_d ) * ddzw(k) &
	4107	) + div * v(k,j,i)
	4108
	4109
	4110	!++
	4111	!-- Statistical Evaluation of v'v'. The factor has to be applied
	4112	!-- for right evaluation when gallilei_trans = .T. .
	4113	! sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
	4114	! + ( flux_n &
	4115	! * ( v_comp - 2.0 * hom(k,1,2,0) ) &
	4116	! / ( v_comp - gv + 1.0E-20 ) &
	4117	! + diss_n &
	4118	! * ABS( v_comp - 2.0 * hom(k,1,2,0) ) &
	4119	! / ( ABS( v_comp - gv ) +1.0E-20 ) ) &
	4120	! * weight_substep(intermediate_timestep_count)
	4121	!
	4122	!-- Statistical Evaluation of w'v'.
	4123	! sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
	4124	! + ( flux_t + diss_t ) &
	4125	! * weight_substep(intermediate_timestep_count)
	4126
	4127	ENDDO
	4128	ENDDO
	4129	ENDDO
	4130	!$acc end kernels
	4131
	4132	!++
	4133	! sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn)
	4134
	4135	END SUBROUTINE advec_v_ws_acc
	4136
	4137
	4138	!------------------------------------------------------------------------------!
[713]	4139	! Advection of w - Call for all grid points
	4140	!------------------------------------------------------------------------------!
	4141	SUBROUTINE advec_w_ws
	4142
	4143	USE arrays_3d
	4144	USE constants
	4145	USE control_parameters
	4146	USE grid_variables
	4147	USE indices
	4148	USE statistics
	4149
	4150	IMPLICIT NONE
	4151
[888]	4152	INTEGER :: i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
	4153	ibit34, ibit35, j, k, k_mm, k_pp, k_ppp, tn = 0
[862]	4154	REAL :: diss_d, div, flux_d, gu, gv, u_comp, v_comp, w_comp
	4155	REAL, DIMENSION(nzb:nzt) :: diss_t, flux_t
	4156	REAL, DIMENSION(nzb+1:nzt) :: diss_n, diss_r, flux_n, flux_r, &
	4157	swap_diss_y_local_w, &
	4158	swap_flux_y_local_w
	4159	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_diss_x_local_w, &
	4160	swap_flux_x_local_w
[713]	4161
	4162	gu = 2.0 * u_gtrans
	4163	gv = 2.0 * v_gtrans
	4164	!
	4165	!-- compute the whole left boundary of the processor domain
	4166	i = nxl
	4167	DO j = nys, nyn
[861]	4168	DO k = nzb+1, nzb_max
	4169
[888]	4170	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
	4171	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
	4172	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
	4173
[861]	4174	u_comp = u(k+1,j,i) + u(k,j,i) - gu
	4175	swap_flux_x_local_w(k,j) = u_comp * ( &
[888]	4176	( 37.0 * ibit29 * adv_mom_5 &
	4177	+ 7.0 * ibit28 * adv_mom_3 &
	4178	+ ibit27 * adv_mom_1 &
	4179	) * &
[861]	4180	( w(k,j,i) + w(k,j,i-1) ) &
[888]	4181	- ( 8.0 * ibit29 * adv_mom_5 &
	4182	+ ibit28 * adv_mom_3 &
	4183	) * &
[861]	4184	( w(k,j,i+1) + w(k,j,i-2) ) &
[888]	4185	+ ( ibit29 * adv_mom_5 &
	4186	) * &
[861]	4187	( w(k,j,i+2) + w(k,j,i-3) ) &
	4188	)
	4189
	4190	swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( &
[888]	4191	( 10.0 * ibit29 * adv_mom_5 &
	4192	+ 3.0 * ibit28 * adv_mom_3 &
	4193	+ ibit27 * adv_mom_1 &
	4194	) * &
[861]	4195	( w(k,j,i) - w(k,j,i-1) ) &
[888]	4196	- ( 5.0 * ibit29 * adv_mom_5 &
	4197	+ ibit28 * adv_mom_3 &
	4198	) * &
[861]	4199	( w(k,j,i+1) - w(k,j,i-2) ) &
[888]	4200	+ ( ibit29 * adv_mom_5 &
	4201	) * &
[861]	4202	( w(k,j,i+2) - w(k,j,i-3) ) &
	4203	)
	4204
	4205	ENDDO
	4206
	4207	DO k = nzb_max+1, nzt
	4208
	4209	u_comp = u(k+1,j,i) + u(k,j,i) - gu
	4210	swap_flux_x_local_w(k,j) = u_comp * ( &
	4211	37.0 * ( w(k,j,i) + w(k,j,i-1) ) &
	4212	- 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) &
	4213	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
	4214	swap_diss_x_local_w(k,j) = - ABS( u_comp ) * ( &
	4215	10.0 * ( w(k,j,i) - w(k,j,i-1) ) &
	4216	- 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) &
	4217	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
	4218
	4219	ENDDO
	4220
[713]	4221	ENDDO
	4222
	4223	DO i = nxl, nxr
	4224
[861]	4225	j = nys
	4226	DO k = nzb+1, nzb_max
[713]	4227
[1221]	4228	ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
[888]	4229	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
	4230	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
	4231
[861]	4232	v_comp = v(k+1,j,i) + v(k,j,i) - gv
	4233	swap_flux_y_local_w(k) = v_comp * ( &
[888]	4234	( 37.0 * ibit32 * adv_mom_5 &
	4235	+ 7.0 * ibit31 * adv_mom_3 &
	4236	+ ibit30 * adv_mom_1 &
	4237	) * &
[861]	4238	( w(k,j,i) + w(k,j-1,i) ) &
[888]	4239	- ( 8.0 * ibit32 * adv_mom_5 &
	4240	+ ibit31 * adv_mom_3 &
	4241	) * &
[861]	4242	( w(k,j+1,i) + w(k,j-2,i) ) &
[888]	4243	+ ( ibit32 * adv_mom_5 &
	4244	) * &
[861]	4245	( w(k,j+2,i) + w(k,j-3,i) ) &
	4246	)
	4247
	4248	swap_diss_y_local_w(k) = - ABS( v_comp ) * ( &
[888]	4249	( 10.0 * ibit32 * adv_mom_5 &
	4250	+ 3.0 * ibit31 * adv_mom_3 &
	4251	+ ibit30 * adv_mom_1 &
	4252	) * &
[861]	4253	( w(k,j,i) - w(k,j-1,i) ) &
[888]	4254	- ( 5.0 * ibit32 * adv_mom_5 &
	4255	+ ibit31 * adv_mom_3 &
	4256	) * &
[861]	4257	( w(k,j+1,i) - w(k,j-2,i) ) &
[888]	4258	+ ( ibit32 * adv_mom_5 &
	4259	) * &
[861]	4260	( w(k,j+2,i) - w(k,j-3,i) ) &
	4261	)
	4262
	4263	ENDDO
	4264
	4265	DO k = nzb_max+1, nzt
	4266
	4267	v_comp = v(k+1,j,i) + v(k,j,i) - gv
	4268	swap_flux_y_local_w(k) = v_comp * ( &
	4269	37.0 * ( w(k,j,i) + w(k,j-1,i) ) &
	4270	- 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) &
	4271	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
	4272	swap_diss_y_local_w(k) = - ABS( v_comp ) * ( &
	4273	10.0 * ( w(k,j,i) - w(k,j-1,i) ) &
	4274	- 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) &
	4275	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
	4276
	4277	ENDDO
	4278
[713]	4279	DO j = nys, nyn
[861]	4280
[862]	4281	!
	4282	!-- The lower flux has to be calculated explicetely for the tendency
	4283	!-- at the first w-level. For topography wall this is done implicitely
	4284	!-- by wall_flags_0.
	4285	k = nzb + 1
	4286	w_comp = w(k,j,i) + w(k-1,j,i)
	4287	flux_t(0) = w_comp * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1
	4288	diss_t(0) = -ABS(w_comp) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1
	4289	flux_d = flux_t(0)
	4290	diss_d = diss_t(0)
[861]	4291
	4292	DO k = nzb+1, nzb_max
	4293
[888]	4294	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
	4295	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
	4296	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
	4297
[861]	4298	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
	4299	flux_r(k) = u_comp * ( &
[888]	4300	( 37.0 * ibit29 * adv_mom_5 &
	4301	+ 7.0 * ibit28 * adv_mom_3 &
	4302	+ ibit27 * adv_mom_1 &
	4303	) * &
[861]	4304	( w(k,j,i+1) + w(k,j,i) ) &
[888]	4305	- ( 8.0 * ibit29 * adv_mom_5 &
	4306	+ ibit28 * adv_mom_3 &
	4307	) * &
[861]	4308	( w(k,j,i+2) + w(k,j,i-1) ) &
[888]	4309	+ ( ibit29 * adv_mom_5 &
	4310	) * &
[861]	4311	( w(k,j,i+3) + w(k,j,i-2) ) &
	4312	)
	4313
	4314	diss_r(k) = - ABS( u_comp ) * ( &
[888]	4315	( 10.0 * ibit29 * adv_mom_5 &
	4316	+ 3.0 * ibit28 * adv_mom_3 &
	4317	+ ibit27 * adv_mom_1 &
	4318	) * &
[861]	4319	( w(k,j,i+1) - w(k,j,i) ) &
[888]	4320	- ( 5.0 * ibit29 * adv_mom_5 &
	4321	+ ibit28 * adv_mom_3 &
	4322	) * &
[861]	4323	( w(k,j,i+2) - w(k,j,i-1) ) &
[888]	4324	+ ( ibit29 * adv_mom_5 &
	4325	) * &
[861]	4326	( w(k,j,i+3) - w(k,j,i-2) ) &
	4327	)
	4328
[1221]	4329	ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
[888]	4330	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
	4331	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
	4332
[861]	4333	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
	4334	flux_n(k) = v_comp * ( &
[888]	4335	( 37.0 * ibit32 * adv_mom_5 &
	4336	+ 7.0 * ibit31 * adv_mom_3 &
	4337	+ ibit30 * adv_mom_1 &
	4338	) * &
[861]	4339	( w(k,j+1,i) + w(k,j,i) ) &
[888]	4340	- ( 8.0 * ibit32 * adv_mom_5 &
	4341	+ ibit31 * adv_mom_3 &
	4342	) * &
[861]	4343	( w(k,j+2,i) + w(k,j-1,i) ) &
[888]	4344	+ ( ibit32 * adv_mom_5 &
	4345	) * &
[861]	4346	( w(k,j+3,i) + w(k,j-2,i) ) &
	4347	)
	4348
	4349	diss_n(k) = - ABS( v_comp ) * ( &
[888]	4350	( 10.0 * ibit32 * adv_mom_5 &
	4351	+ 3.0 * ibit31 * adv_mom_3 &
	4352	+ ibit30 * adv_mom_1 &
	4353	) * &
[861]	4354	( w(k,j+1,i) - w(k,j,i) ) &
[888]	4355	- ( 5.0 * ibit32 * adv_mom_5 &
	4356	+ ibit31 * adv_mom_3 &
	4357	) * &
[861]	4358	( w(k,j+2,i) - w(k,j-1,i) ) &
[888]	4359	+ ( ibit32 * adv_mom_5 &
	4360	) * &
[861]	4361	( w(k,j+3,i) - w(k,j-2,i) ) &
	4362	)
[713]	4363	!
[861]	4364	!-- k index has to be modified near bottom and top, else array
	4365	!-- subscripts will be exceeded.
[1221]	4366	ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
	4367	ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
	4368	ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
[713]	4369
[888]	4370	k_ppp = k + 3 * ibit35
	4371	k_pp = k + 2 * ( 1 - ibit33 )
	4372	k_mm = k - 2 * ibit35
	4373
[861]	4374	w_comp = w(k+1,j,i) + w(k,j,i)
	4375	flux_t(k) = w_comp * ( &
[888]	4376	( 37.0 * ibit35 * adv_mom_5 &
	4377	+ 7.0 * ibit34 * adv_mom_3 &
	4378	+ ibit33 * adv_mom_1 &
	4379	) * &
[861]	4380	( w(k+1,j,i) + w(k,j,i) ) &
[888]	4381	- ( 8.0 * ibit35 * adv_mom_5 &
	4382	+ ibit34 * adv_mom_3 &
	4383	) * &
[861]	4384	( w(k_pp,j,i) + w(k-1,j,i) ) &
[888]	4385	+ ( ibit35 * adv_mom_5 &
	4386	) * &
[861]	4387	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
	4388	)
	4389
	4390	diss_t(k) = - ABS( w_comp ) * ( &
[888]	4391	( 10.0 * ibit35 * adv_mom_5 &
	4392	+ 3.0 * ibit34 * adv_mom_3 &
	4393	+ ibit33 * adv_mom_1 &
	4394	) * &
[861]	4395	( w(k+1,j,i) - w(k,j,i) ) &
[888]	4396	- ( 5.0 * ibit35 * adv_mom_5 &
	4397	+ ibit34 * adv_mom_3 &
	4398	) * &
[861]	4399	( w(k_pp,j,i) - w(k-1,j,i) ) &
[888]	4400	+ ( ibit35 * adv_mom_5 &
	4401	) * &
[861]	4402	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
	4403	)
[713]	4404	!
[861]	4405	!-- Calculate the divergence of the velocity field. A respective
	4406	!-- correction is needed to overcome numerical instabilities caused
	4407	!-- by a not sufficient reduction of divergences near topography.
	4408	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
	4409	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
	4410	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) &
	4411	* ddzu(k+1) &
	4412	) * 0.5
[713]	4413
[861]	4414	tend(k,j,i) = tend(k,j,i) - ( &
	4415	( flux_r(k) + diss_r(k) &
	4416	- swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) &
	4417	) * ddx &
	4418	+ ( flux_n(k) + diss_n(k) &
	4419	- swap_flux_y_local_w(k) - swap_diss_y_local_w(k) &
	4420	) * ddy &
	4421	+ ( flux_t(k) + diss_t(k) &
	4422	- flux_d - diss_d &
	4423	) * ddzu(k+1) &
	4424	) + div * w(k,j,i)
[713]	4425
[861]	4426	swap_flux_x_local_w(k,j) = flux_r(k)
	4427	swap_diss_x_local_w(k,j) = diss_r(k)
	4428	swap_flux_y_local_w(k) = flux_n(k)
	4429	swap_diss_y_local_w(k) = diss_n(k)
	4430	flux_d = flux_t(k)
	4431	diss_d = diss_t(k)
	4432
	4433	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
	4434	+ ( flux_t(k) + diss_t(k) ) &
	4435	* weight_substep(intermediate_timestep_count)
	4436
	4437	ENDDO
	4438
	4439	DO k = nzb_max+1, nzt
	4440
	4441	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
	4442	flux_r(k) = u_comp * ( &
	4443	37.0 * ( w(k,j,i+1) + w(k,j,i) ) &
	4444	- 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) &
	4445	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
	4446
	4447	diss_r(k) = - ABS( u_comp ) * ( &
	4448	10.0 * ( w(k,j,i+1) - w(k,j,i) ) &
	4449	- 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) &
	4450	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
	4451
	4452	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
	4453	flux_n(k) = v_comp * ( &
	4454	37.0 * ( w(k,j+1,i) + w(k,j,i) ) &
	4455	- 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) &
	4456	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
	4457
	4458	diss_n(k) = - ABS( v_comp ) * ( &
	4459	10.0 * ( w(k,j+1,i) - w(k,j,i) ) &
	4460	- 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) &
	4461	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
[713]	4462	!
[861]	4463	!-- k index has to be modified near bottom and top, else array
	4464	!-- subscripts will be exceeded.
[1221]	4465	ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
	4466	ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
	4467	ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
[713]	4468
[888]	4469	k_ppp = k + 3 * ibit35
	4470	k_pp = k + 2 * ( 1 - ibit33 )
	4471	k_mm = k - 2 * ibit35
	4472
[861]	4473	w_comp = w(k+1,j,i) + w(k,j,i)
	4474	flux_t(k) = w_comp * ( &
[888]	4475	( 37.0 * ibit35 * adv_mom_5 &
	4476	+ 7.0 * ibit34 * adv_mom_3 &
	4477	+ ibit33 * adv_mom_1 &
	4478	) * &
[861]	4479	( w(k+1,j,i) + w(k,j,i) ) &
[888]	4480	- ( 8.0 * ibit35 * adv_mom_5 &
	4481	+ ibit34 * adv_mom_3 &
	4482	) * &
[861]	4483	( w(k_pp,j,i) + w(k-1,j,i) ) &
[888]	4484	+ ( ibit35 * adv_mom_5 &
	4485	) * &
[861]	4486	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
[888]	4487	)
[861]	4488
	4489	diss_t(k) = - ABS( w_comp ) * ( &
[888]	4490	( 10.0 * ibit35 * adv_mom_5 &
	4491	+ 3.0 * ibit34 * adv_mom_3 &
	4492	+ ibit33 * adv_mom_1 &
	4493	) * &
[861]	4494	( w(k+1,j,i) - w(k,j,i) ) &
[888]	4495	- ( 5.0 * ibit35 * adv_mom_5 &
	4496	+ ibit34 * adv_mom_3 &
	4497	) * &
[861]	4498	( w(k_pp,j,i) - w(k-1,j,i) ) &
[888]	4499	+ ( ibit35 * adv_mom_5 &
	4500	) * &
[861]	4501	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
[888]	4502	)
[713]	4503	!
[861]	4504	!-- Calculate the divergence of the velocity field. A respective
	4505	!-- correction is needed to overcome numerical instabilities caused
	4506	!-- by a not sufficient reduction of divergences near topography.
	4507	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
	4508	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
	4509	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) &
	4510	* ddzu(k+1) &
	4511	) * 0.5
	4512
	4513	tend(k,j,i) = tend(k,j,i) - ( &
	4514	( flux_r(k) + diss_r(k) &
	4515	- swap_flux_x_local_w(k,j) - swap_diss_x_local_w(k,j) &
	4516	) * ddx &
	4517	+ ( flux_n(k) + diss_n(k) &
	4518	- swap_flux_y_local_w(k) - swap_diss_y_local_w(k) &
	4519	) * ddy &
	4520	+ ( flux_t(k) + diss_t(k) &
	4521	- flux_d - diss_d &
	4522	) * ddzu(k+1) &
	4523	) + div * w(k,j,i)
	4524
	4525	swap_flux_x_local_w(k,j) = flux_r(k)
	4526	swap_diss_x_local_w(k,j) = diss_r(k)
	4527	swap_flux_y_local_w(k) = flux_n(k)
	4528	swap_diss_y_local_w(k) = diss_n(k)
	4529	flux_d = flux_t(k)
	4530	diss_d = diss_t(k)
	4531
	4532	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
	4533	+ ( flux_t(k) + diss_t(k) ) &
[743]	4534	* weight_substep(intermediate_timestep_count)
[713]	4535
[861]	4536	ENDDO
[713]	4537	ENDDO
	4538	ENDDO
	4539
	4540	END SUBROUTINE advec_w_ws
	4541
[1015]	4542
	4543	!------------------------------------------------------------------------------!
	4544	! Advection of w - Call for all grid points - accelerator version
	4545	!------------------------------------------------------------------------------!
	4546	SUBROUTINE advec_w_ws_acc
	4547
	4548	USE arrays_3d
	4549	USE constants
	4550	USE control_parameters
	4551	USE grid_variables
	4552	USE indices
	4553	USE statistics
	4554
	4555	IMPLICIT NONE
	4556
	4557	INTEGER :: i, ibit27, ibit28, ibit29, ibit30, ibit31, ibit32, ibit33, &
	4558	ibit34, ibit35, j, k, k_mmm, k_mm, k_pp, k_ppp, tn = 0
	4559
	4560	REAL :: diss_d, diss_l, diss_n, diss_r, diss_s, diss_t, div, &
	4561	flux_d, flux_l, flux_n, flux_r, flux_s, flux_t, gu, gv, &
	4562	u_comp, u_comp_l, v_comp, v_comp_s, w_comp
	4563
	4564	gu = 2.0 * u_gtrans
	4565	gv = 2.0 * v_gtrans
	4566
[1221]	4567
[1015]	4568	!
	4569	!-- Computation of fluxes and tendency terms
[1221]	4570	!$acc kernels present( ddzu, tend, u, v, w, wall_flags_0, wall_flags_00 )
[1015]	4571	!$acc loop
[1128]	4572	DO i = i_left, i_right
	4573	DO j = j_south, j_north
[1015]	4574	!$acc loop vector( 32 )
	4575	DO k = nzb+1, nzt
	4576
[1221]	4577	ibit27 = IBITS(wall_flags_0(k,j,i),27,1)
	4578	ibit28 = IBITS(wall_flags_0(k,j,i),28,1)
[1015]	4579	ibit29 = IBITS(wall_flags_0(k,j,i),29,1)
	4580
	4581	u_comp_l = u(k+1,j,i) + u(k,j,i) - gu
	4582	flux_l = u_comp_l * ( &
	4583	( 37.0 * ibit29 * adv_mom_5 &
	4584	+ 7.0 * ibit28 * adv_mom_3 &
	4585	+ ibit27 * adv_mom_1 &
	4586	) * &
	4587	( w(k,j,i) + w(k,j,i-1) ) &
	4588	- ( 8.0 * ibit29 * adv_mom_5 &
	4589	+ ibit28 * adv_mom_3 &
	4590	) * &
	4591	( w(k,j,i+1) + w(k,j,i-2) ) &
	4592	+ ( ibit29 * adv_mom_5 &
	4593	) * &
	4594	( w(k,j,i+2) + w(k,j,i-3) ) &
	4595	)
	4596
	4597	diss_l = - ABS( u_comp_l ) * ( &
	4598	( 10.0 * ibit29 * adv_mom_5 &
	4599	+ 3.0 * ibit28 * adv_mom_3 &
	4600	+ ibit27 * adv_mom_1 &
	4601	) * &
	4602	( w(k,j,i) - w(k,j,i-1) ) &
	4603	- ( 5.0 * ibit29 * adv_mom_5 &
	4604	+ ibit28 * adv_mom_3 &
	4605	) * &
	4606	( w(k,j,i+1) - w(k,j,i-2) ) &
	4607	+ ( ibit29 * adv_mom_5 &
	4608	) * &
	4609	( w(k,j,i+2) - w(k,j,i-3) ) &
	4610	)
	4611
	4612	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
	4613	flux_r = u_comp * ( &
	4614	( 37.0 * ibit29 * adv_mom_5 &
	4615	+ 7.0 * ibit28 * adv_mom_3 &
	4616	+ ibit27 * adv_mom_1 &
	4617	) * &
	4618	( w(k,j,i+1) + w(k,j,i) ) &
	4619	- ( 8.0 * ibit29 * adv_mom_5 &
	4620	+ ibit28 * adv_mom_3 &
	4621	) * &
	4622	( w(k,j,i+2) + w(k,j,i-1) ) &
	4623	+ ( ibit29 * adv_mom_5 &
	4624	) * &
	4625	( w(k,j,i+3) + w(k,j,i-2) ) &
	4626	)
	4627
	4628	diss_r = - ABS( u_comp ) * ( &
	4629	( 10.0 * ibit29 * adv_mom_5 &
	4630	+ 3.0 * ibit28 * adv_mom_3 &
	4631	+ ibit27 * adv_mom_1 &
	4632	) * &
	4633	( w(k,j,i+1) - w(k,j,i) ) &
	4634	- ( 5.0 * ibit29 * adv_mom_5 &
	4635	+ ibit28 * adv_mom_3 &
	4636	) * &
	4637	( w(k,j,i+2) - w(k,j,i-1) ) &
	4638	+ ( ibit29 * adv_mom_5 &
	4639	) * &
	4640	( w(k,j,i+3) - w(k,j,i-2) ) &
	4641	)
[1221]	4642	ibit32 = IBITS(wall_flags_00(k,j,i),0,1)
[1015]	4643	ibit31 = IBITS(wall_flags_0(k,j,i),31,1)
	4644	ibit30 = IBITS(wall_flags_0(k,j,i),30,1)
	4645
	4646	v_comp_s = v(k+1,j,i) + v(k,j,i) - gv
	4647	flux_s = v_comp_s * ( &
	4648	( 37.0 * ibit32 * adv_mom_5 &
	4649	+ 7.0 * ibit31 * adv_mom_3 &
	4650	+ ibit30 * adv_mom_1 &
	4651	) * &
	4652	( w(k,j,i) + w(k,j-1,i) ) &
	4653	- ( 8.0 * ibit32 * adv_mom_5 &
	4654	+ ibit31 * adv_mom_3 &
	4655	) * &
	4656	( w(k,j+1,i) + w(k,j-2,i) ) &
	4657	+ ( ibit32 * adv_mom_5 &
	4658	) * &
	4659	( w(k,j+2,i) + w(k,j-3,i) ) &
	4660	)
	4661
	4662	diss_s = - ABS( v_comp_s ) * ( &
	4663	( 10.0 * ibit32 * adv_mom_5 &
	4664	+ 3.0 * ibit31 * adv_mom_3 &
	4665	+ ibit30 * adv_mom_1 &
	4666	) * &
	4667	( w(k,j,i) - w(k,j-1,i) ) &
	4668	- ( 5.0 * ibit32 * adv_mom_5 &
	4669	+ ibit31 * adv_mom_3 &
	4670	) * &
	4671	( w(k,j+1,i) - w(k,j-2,i) ) &
	4672	+ ( ibit32 * adv_mom_5 &
	4673	) * &
	4674	( w(k,j+2,i) - w(k,j-3,i) ) &
	4675	)
	4676
	4677	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
	4678	flux_n = v_comp * ( &
	4679	( 37.0 * ibit32 * adv_mom_5 &
	4680	+ 7.0 * ibit31 * adv_mom_3 &
	4681	+ ibit30 * adv_mom_1 &
	4682	) * &
	4683	( w(k,j+1,i) + w(k,j,i) ) &
	4684	- ( 8.0 * ibit32 * adv_mom_5 &
	4685	+ ibit31 * adv_mom_3 &
	4686	) * &
	4687	( w(k,j+2,i) + w(k,j-1,i) ) &
	4688	+ ( ibit32 * adv_mom_5 &
	4689	) * &
	4690	( w(k,j+3,i) + w(k,j-2,i) ) &
	4691	)
	4692
	4693	diss_n = - ABS( v_comp ) * ( &
	4694	( 10.0 * ibit32 * adv_mom_5 &
	4695	+ 3.0 * ibit31 * adv_mom_3 &
	4696	+ ibit30 * adv_mom_1 &
	4697	) * &
	4698	( w(k,j+1,i) - w(k,j,i) ) &
	4699	- ( 5.0 * ibit32 * adv_mom_5 &
	4700	+ ibit31 * adv_mom_3 &
	4701	) * &
	4702	( w(k,j+2,i) - w(k,j-1,i) ) &
	4703	+ ( ibit32 * adv_mom_5 &
	4704	) * &
	4705	( w(k,j+3,i) - w(k,j-2,i) ) &
	4706	)
	4707
[1221]	4708	ibit35 = IBITS(wall_flags_00(k-1,j,i),3,1)
	4709	ibit34 = IBITS(wall_flags_00(k-1,j,i),2,1)
	4710	ibit33 = IBITS(wall_flags_00(k-1,j,i),1,1)
[1015]	4711
[1027]	4712	k_pp = k + 2 * ibit35
	4713	k_mm = k - 2 * ( ibit34 + ibit35 )
[1015]	4714	k_mmm = k - 3 * ibit35
	4715
	4716	w_comp = w(k,j,i) + w(k-1,j,i)
	4717	flux_d = w_comp * ( &
	4718	( 37.0 * ibit35 * adv_mom_5 &
	4719	+ 7.0 * ibit34 * adv_mom_3 &
	4720	+ ibit33 * adv_mom_1 &
	4721	) * &
	4722	( w(k,j,i) + w(k-1,j,i) ) &
	4723	- ( 8.0 * ibit35 * adv_mom_5 &
	4724	+ ibit34 * adv_mom_3 &
	4725	) * &
	4726	( w(k+1,j,i) + w(k_mm,j,i) ) &
	4727	+ ( ibit35 * adv_mom_5 &
	4728	) * &
	4729	( w(k_pp,j,i) + w(k_mmm,j,i) ) &
	4730	)
	4731
	4732	diss_d = - ABS( w_comp ) * ( &
	4733	( 10.0 * ibit35 * adv_mom_5 &
	4734	+ 3.0 * ibit34 * adv_mom_3 &
	4735	+ ibit33 * adv_mom_1 &
	4736	) * &
	4737	( w(k,j,i) - w(k-1,j,i) ) &
	4738	- ( 5.0 * ibit35 * adv_mom_5 &
	4739	+ ibit34 * adv_mom_3 &
	4740	) * &
	4741	( w(k+1,j,i) - w(k_mm,j,i) ) &
	4742	+ ( ibit35 * adv_mom_5 &
	4743	) * &
	4744	( w(k_pp,j,i) - w(k_mmm,j,i) ) &
	4745	)
	4746
	4747	!
	4748	!-- k index has to be modified near bottom and top, else array
	4749	!-- subscripts will be exceeded.
[1221]	4750	ibit35 = IBITS(wall_flags_00(k,j,i),3,1)
	4751	ibit34 = IBITS(wall_flags_00(k,j,i),2,1)
	4752	ibit33 = IBITS(wall_flags_00(k,j,i),1,1)
[1015]	4753
	4754	k_ppp = k + 3 * ibit35
	4755	k_pp = k + 2 * ( 1 - ibit33 )
	4756	k_mm = k - 2 * ibit35
	4757
	4758	w_comp = w(k+1,j,i) + w(k,j,i)
	4759	flux_t = w_comp * ( &
	4760	( 37.0 * ibit35 * adv_mom_5 &
	4761	+ 7.0 * ibit34 * adv_mom_3 &
	4762	+ ibit33 * adv_mom_1 &
	4763	) * &
	4764	( w(k+1,j,i) + w(k,j,i) ) &
	4765	- ( 8.0 * ibit35 * adv_mom_5 &
	4766	+ ibit34 * adv_mom_3 &
	4767	) * &
	4768	( w(k_pp,j,i) + w(k-1,j,i) ) &
	4769	+ ( ibit35 * adv_mom_5 &
	4770	) * &
	4771	( w(k_ppp,j,i) + w(k_mm,j,i) ) &
	4772	)
	4773
	4774	diss_t = - ABS( w_comp ) * ( &
	4775	( 10.0 * ibit35 * adv_mom_5 &
	4776	+ 3.0 * ibit34 * adv_mom_3 &
	4777	+ ibit33 * adv_mom_1 &
	4778	) * &
	4779	( w(k+1,j,i) - w(k,j,i) ) &
	4780	- ( 5.0 * ibit35 * adv_mom_5 &
	4781	+ ibit34 * adv_mom_3 &
	4782	) * &
	4783	( w(k_pp,j,i) - w(k-1,j,i) ) &
	4784	+ ( ibit35 * adv_mom_5 &
	4785	) * &
	4786	( w(k_ppp,j,i) - w(k_mm,j,i) ) &
	4787	)
	4788	!
	4789	!-- Calculate the divergence of the velocity field. A respective
	4790	!-- correction is needed to overcome numerical instabilities caused
	4791	!-- by a not sufficient reduction of divergences near topography.
	4792	div = ( ( u_comp + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &
	4793	+ ( v_comp + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &
	4794	+ ( w_comp - ( w(k,j,i) + w(k-1,j,i) ) ) &
	4795	* ddzu(k+1) &
	4796	) * 0.5
	4797
	4798	tend(k,j,i) = - ( &
	4799	( flux_r + diss_r - flux_l - diss_l ) * ddx &
	4800	+ ( flux_n + diss_n - flux_s - diss_s ) * ddy &
	4801	+ ( flux_t + diss_t - flux_d - diss_d ) * ddzu(k+1) &
	4802	) + div * w(k,j,i)
	4803
	4804
	4805	!++
	4806	!-- Statistical Evaluation of w'w'.
	4807	! sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
	4808	! + ( flux_t + diss_t ) &
	4809	! * weight_substep(intermediate_timestep_count)
	4810
	4811	ENDDO
	4812	ENDDO
	4813	ENDDO
	4814	!$acc end kernels
	4815
	4816	END SUBROUTINE advec_w_ws_acc
	4817
[713]	4818	END MODULE advec_ws

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

Context Navigation

source: palm/trunk/SOURCE/advec_ws.f90 @ 1221

Download in other formats: