Home

Context Navigation

advec_ws.f90 @ 857

Last change on this file since 857 was 857, checked in by suehring, 12 years ago
last commit documented
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 160.6 KB

Line
1	MODULE advec_ws
2
3	!-----------------------------------------------------------------------------!
4	! Current revisions:
5	! ------------------
6	!
7	! Former revisions:
8	! -----------------
9	! $Id: advec_ws.f90 857 2012-03-20 18:30:38Z suehring $
10	!
11	! 856 2012-03-20 18:29:17Z suehring
12	! Bug concerning numerical dissipation at the first grid level. Function
13	! diss_2nd always returns zero. Because diss_2nd will be replaced by numerical
14	! dissipation which is more consistent to the WS-schemes, the return value
15	! is set to zero for now.
16	!
17	! 801 2012-01-10 17:30:36Z suehring
18	! Bugfix concerning OpenMP parallelization. Summation of sums_wsus_ws_l,
19	! sums_wsvs_ws_l, sums_us2_ws_l, sums_vs2_ws_l, sums_ws2_ws_l, sums_wspts_ws_l,
20	! sums_wsqs_ws_l, sums_wssas_ws_l is now thread-safe by adding an additional
21	! dimension.
22	!
23	! 743 2011-08-18 16:10:16Z suehring
24	! Evaluation of turbulent fluxes with WS-scheme only for the whole model
25	! domain. Therefor dimension of arrays needed for statistical evaluation
26	! decreased.
27	!
28	! 736 2011-08-17 14:13:26Z suehring
29	! Bugfix concerning OpenMP parallelization. i_omp introduced, because first
30	! index where fluxes on left side have to be calculated explicitly is
31	! different on each thread. Furthermore the swapping of fluxes is now
32	! thread-safe by adding an additional dimension.
33	!
34	! 713 2011-03-30 14:21:21Z suehring
35	! File reformatted.
36	! Bugfix in vertical advection of w concerning the optimized version for
37	! vector architecture.
38	! Constants adv_mom_3, adv_mom_5, adv_sca_5, adv_sca_3 reformulated as
39	! broken numbers.
40	!
41	! 709 2011-03-30 09:31:40Z raasch
42	! formatting adjustments
43	!
44	! 705 2011-03-25 11:21:43 Z suehring $
45	! Bugfix in declaration of logicals concerning outflow boundaries.
46	!
47	! 411 2009-12-11 12:31:43 Z suehring
48	! Allocation of weight_substep moved to init_3d_model.
49	! Declaration of ws_scheme_sca and ws_scheme_mom moved to check_parameters.
50	! Setting bc for the horizontal velocity variances added (moved from
51	! flow_statistics).
52	!
53	! Initial revision
54	!
55	! 411 2009-12-11 12:31:43 Z suehring
56	!
57	! Description:
58	! ------------
59	! Advection scheme for scalars and momentum using the flux formulation of
60	! Wicker and Skamarock 5th order. Additionally the module contains of a
61	! routine using for initialisation and steering of the statical evaluation.
62	! The computation of turbulent fluxes takes place inside the advection
63	! routines.
64	! In case of vector architectures Dirichlet and Radiation boundary conditions
65	! are outstanding and not available.
66	! A further routine local_diss_ij is available (next weeks) and is used if a
67	! control of dissipative fluxes is desired.
68	!
69	! OUTSTANDING: - Dirichlet and Radiation boundary conditions for
70	! vector architectures
71	! - dissipation control for cache architectures ( next weeks )
72	! - Topography ( next weeks )
73	!-----------------------------------------------------------------------------!
74
75	PRIVATE
76	PUBLIC advec_s_ws, advec_u_ws, advec_v_ws, advec_w_ws, &
77	local_diss, ws_init, ws_statistics
78
79	INTERFACE ws_init
80	MODULE PROCEDURE ws_init
81	END INTERFACE ws_init
82
83	INTERFACE ws_statistics
84	MODULE PROCEDURE ws_statistics
85	END INTERFACE ws_statistics
86
87	INTERFACE advec_s_ws
88	MODULE PROCEDURE advec_s_ws
89	MODULE PROCEDURE advec_s_ws_ij
90	END INTERFACE advec_s_ws
91
92	INTERFACE advec_u_ws
93	MODULE PROCEDURE advec_u_ws
94	MODULE PROCEDURE advec_u_ws_ij
95	END INTERFACE advec_u_ws
96
97	INTERFACE advec_v_ws
98	MODULE PROCEDURE advec_v_ws
99	MODULE PROCEDURE advec_v_ws_ij
100	END INTERFACE advec_v_ws
101
102	INTERFACE advec_w_ws
103	MODULE PROCEDURE advec_w_ws
104	MODULE PROCEDURE advec_w_ws_ij
105	END INTERFACE advec_w_ws
106
107	INTERFACE local_diss
108	MODULE PROCEDURE local_diss
109	MODULE PROCEDURE local_diss_ij
110	END INTERFACE local_diss
111
112	CONTAINS
113
114
115	!------------------------------------------------------------------------------!
116	! Initialization of WS-scheme
117	!------------------------------------------------------------------------------!
118	SUBROUTINE ws_init
119
120	USE arrays_3d
121	USE constants
122	USE control_parameters
123	USE indices
124	USE pegrid
125	USE statistics
126
127	!
128	!-- Allocate arrays needed for dissipation control.
129	IF ( dissipation_control ) THEN
130	! ALLOCATE(var_x(nzb+1:nzt,nys:nyn,nxl:nxr), &
131	! var_y(nzb+1:nzt,nys:nyn,nxl:nxr), &
132	! var_z(nzb+1:nzt,nys:nyn,nxl:nxr), &
133	! gamma_x(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
134	! gamma_y(nzb:nzt+1,nysg:nyng,nxlg:nxrg), &
135	! gamma_z(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
136	ENDIF
137
138	!
139	!-- Set the appropriate factors for scalar and momentum advection.
140	adv_sca_5 = 1./60.
141	adv_sca_3 = 1./12.
142	adv_mom_5 = 1./120.
143	adv_mom_3 = 1./24.
144
145	!
146	!-- Arrays needed for statical evaluation of fluxes.
147	IF ( ws_scheme_mom ) THEN
148
149	ALLOCATE( sums_wsus_ws_l(nzb:nzt+1,0:threads_per_task-1), &
150	sums_wsvs_ws_l(nzb:nzt+1,0:threads_per_task-1), &
151	sums_us2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
152	sums_vs2_ws_l(nzb:nzt+1,0:threads_per_task-1), &
153	sums_ws2_ws_l(nzb:nzt+1,0:threads_per_task-1) )
154
155	sums_wsus_ws_l = 0.0
156	sums_wsvs_ws_l = 0.0
157	sums_us2_ws_l = 0.0
158	sums_vs2_ws_l = 0.0
159	sums_ws2_ws_l = 0.0
160
161	ENDIF
162
163	IF ( ws_scheme_sca ) THEN
164
165	ALLOCATE( sums_wspts_ws_l(nzb:nzt+1,0:threads_per_task-1) )
166	sums_wspts_ws_l = 0.0
167
168	IF ( humidity .OR. passive_scalar ) THEN
169	ALLOCATE( sums_wsqs_ws_l(nzb:nzt+1,0:threads_per_task-1) )
170	sums_wsqs_ws_l = 0.0
171	ENDIF
172
173	IF ( ocean ) THEN
174	ALLOCATE( sums_wssas_ws_l(nzb:nzt+1,0:threads_per_task-1) )
175	sums_wssas_ws_l = 0.0
176	ENDIF
177
178	ENDIF
179
180	!
181	!-- Arrays needed for reasons of speed optimization for cache and noopt
182	!-- version. For the vector version the buffer arrays are not necessary,
183	!-- because the the fluxes can swapped directly inside the loops of the
184	!-- advection routines.
185	IF ( loop_optimization /= 'vector' ) THEN
186
187	IF ( ws_scheme_mom ) THEN
188
189	ALLOCATE( flux_s_u(nzb+1:nzt,0:threads_per_task-1), &
190	flux_s_v(nzb+1:nzt,0:threads_per_task-1), &
191	flux_s_w(nzb+1:nzt,0:threads_per_task-1), &
192	diss_s_u(nzb+1:nzt,0:threads_per_task-1), &
193	diss_s_v(nzb+1:nzt,0:threads_per_task-1), &
194	diss_s_w(nzb+1:nzt,0:threads_per_task-1) )
195	ALLOCATE( flux_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
196	flux_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
197	flux_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
198	diss_l_u(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
199	diss_l_v(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
200	diss_l_w(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
201
202	ENDIF
203
204	IF ( ws_scheme_sca ) THEN
205
206	ALLOCATE( flux_s_pt(nzb+1:nzt,0:threads_per_task-1), &
207	flux_s_e(nzb+1:nzt,0:threads_per_task-1), &
208	diss_s_pt(nzb+1:nzt,0:threads_per_task-1), &
209	diss_s_e(nzb+1:nzt,0:threads_per_task-1) )
210	ALLOCATE( flux_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
211	flux_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
212	diss_l_pt(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
213	diss_l_e(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
214
215	IF ( humidity .OR. passive_scalar ) THEN
216	ALLOCATE( flux_s_q(nzb+1:nzt,0:threads_per_task-1), &
217	diss_s_q(nzb+1:nzt,0:threads_per_task-1) )
218	ALLOCATE( flux_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
219	diss_l_q(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
220	ENDIF
221
222	IF ( ocean ) THEN
223	ALLOCATE( flux_s_sa(nzb+1:nzt,0:threads_per_task-1), &
224	diss_s_sa(nzb+1:nzt,0:threads_per_task-1) )
225	ALLOCATE( flux_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1), &
226	diss_l_sa(nzb+1:nzt,nys:nyn,0:threads_per_task-1) )
227	ENDIF
228
229	ENDIF
230
231	ENDIF
232
233	!
234	!-- Determine the flags where the order of the scheme for horizontal
235	!-- advection has to be degraded.
236	ALLOCATE( boundary_flags(nys:nyn,nxl:nxr) )
237	DO i = nxl, nxr
238	DO j = nys, nyn
239
240	boundary_flags(j,i) = 0
241	IF ( outflow_l ) THEN
242	IF ( i == nxlu ) boundary_flags(j,i) = 5
243	IF ( i == nxl ) boundary_flags(j,i) = 6
244	ELSEIF ( outflow_r ) THEN
245	IF ( i == nxr-1 ) boundary_flags(j,i) = 1
246	IF ( i == nxr ) boundary_flags(j,i) = 2
247	ELSEIF ( outflow_n ) THEN
248	IF ( j == nyn-1 ) boundary_flags(j,i) = 3
249	IF ( j == nyn ) boundary_flags(j,i) = 4
250	ELSEIF ( outflow_s ) THEN
251	IF ( j == nysv ) boundary_flags(j,i) = 7
252	IF ( j == nys ) boundary_flags(j,i) = 8
253	ENDIF
254
255	ENDDO
256	ENDDO
257
258	END SUBROUTINE ws_init
259
260
261	!------------------------------------------------------------------------------!
262	! Initialize variables used for storing statistic qauntities (fluxes, variances)
263	!------------------------------------------------------------------------------!
264	SUBROUTINE ws_statistics
265
266	USE control_parameters
267	USE statistics
268
269	IMPLICIT NONE
270
271	!
272	!-- The arrays needed for statistical evaluation are set to to 0 at the
273	!-- begin of prognostic_equations.
274	IF ( ws_scheme_mom ) THEN
275	sums_wsus_ws_l = 0.0
276	sums_wsvs_ws_l = 0.0
277	sums_us2_ws_l = 0.0
278	sums_vs2_ws_l = 0.0
279	sums_ws2_ws_l = 0.0
280	ENDIF
281
282	IF ( ws_scheme_sca ) THEN
283	sums_wspts_ws_l = 0.0
284	IF ( humidity .OR. passive_scalar ) sums_wsqs_ws_l = 0.0
285	IF ( ocean ) sums_wssas_ws_l = 0.0
286
287	ENDIF
288
289	END SUBROUTINE ws_statistics
290
291
292	!------------------------------------------------------------------------------!
293	! Scalar advection - Call for grid point i,j
294	!------------------------------------------------------------------------------!
295	SUBROUTINE advec_s_ws_ij( i, j, sk, sk_char,swap_flux_y_local, &
296	swap_diss_y_local, swap_flux_x_local, &
297	swap_diss_x_local, i_omp, tn )
298
299	USE arrays_3d
300	USE constants
301	USE control_parameters
302	USE grid_variables
303	USE indices
304	USE pegrid
305	USE statistics
306
307	IMPLICIT NONE
308
309	INTEGER :: i, i_omp, j, k, tn
310	LOGICAL :: degraded_l, degraded_s
311	REAL :: flux_d, diss_d, u_comp, v_comp
312	REAL, DIMENSION(:,:,:), POINTER :: sk
313	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, &
314	flux_n, diss_n
315	REAL, DIMENSION(nzb+1:nzt,0:threads_per_task-1) :: swap_flux_y_local, &
316	swap_diss_y_local
317	REAL, DIMENSION(nzb+1:nzt,nys:nyn,0:threads_per_task-1) :: &
318	swap_flux_x_local, &
319	swap_diss_x_local
320	CHARACTER (LEN = *), INTENT(IN) :: sk_char
321
322
323	degraded_l = .FALSE.
324	degraded_s = .FALSE.
325
326	IF ( boundary_flags(j,i) /= 0 ) THEN
327	!
328	!-- Degrade the order for Dirichlet bc. at the outflow boundary
329	SELECT CASE ( boundary_flags(j,i) )
330
331	CASE ( 1 )
332
333	DO k = nzb_s_inner(j,i)+1, nzt
334	u_comp = u(k,j,i+1) - u_gtrans
335	flux_r(k) = u_comp * ( &
336	7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
337	- ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
338	diss_r(k) = -ABS( u_comp ) * ( &
339	3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
340	- ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
341	ENDDO
342
343	CASE ( 2 )
344
345	DO k = nzb_s_inner(j,i)+1, nzt
346	u_comp = u(k,j,i+1) - u_gtrans
347	flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
348	diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), sk(k,j,i), &
349	sk(k,j,i-1), sk(k,j,i-2), u_comp, &
350	0.5, ddx )
351	ENDDO
352
353	CASE ( 3 )
354
355	DO k = nzb_s_inner(j,i)+1, nzt
356	v_comp = v(k,j+1,i) - v_gtrans
357	flux_n(k) = v_comp * ( &
358	7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
359	- ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
360	diss_n(k) = -ABS( v_comp ) * ( &
361	3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
362	- ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
363	ENDDO
364
365	CASE ( 4 )
366
367	DO k = nzb_s_inner(j,i)+1, nzt
368	v_comp = v(k,j+1,i) - v_gtrans
369	flux_n(k) = v_comp* ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
370	diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), sk(k,j,i), &
371	sk(k,j-1,i), sk(k,j-2,i), v_comp, &
372	0.5, ddy )
373	ENDDO
374
375	CASE ( 5 )
376
377	DO k = nzb_w_inner(j,i)+1, nzt
378	u_comp = u(k,j,i+1) - u_gtrans
379	flux_r(k) = u_comp * ( &
380	7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
381	- ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
382	diss_r(k) = -ABS( u_comp ) * ( &
383	3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
384	- ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
385	ENDDO
386
387	CASE ( 6 )
388
389	DO k = nzb_s_inner(j,i)+1, nzt
390	u_comp = u(k,j,i+1) - u_gtrans
391	flux_r(k) = u_comp * ( &
392	7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
393	- ( sk(k,j,i+2) + sk(k,j,i-1) ) ) * adv_sca_3
394	diss_r(k) = -ABS( u_comp ) * ( &
395	3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
396	- ( sk(k,j,i+2) - sk(k,j,i-1) ) ) * adv_sca_3
397	!
398	!-- Compute leftside fluxes for the left boundary of PE domain
399	u_comp = u(k,j,i) - u_gtrans
400	swap_flux_x_local(k,j,tn) = u_comp * ( &
401	sk(k,j,i) + sk(k,j,i-1) ) * 0.5
402	swap_diss_x_local(k,j,tn) = diss_2nd( sk(k,j,i+2),sk(k,j,i+1), &
403	sk(k,j,i), sk(k,j,i-1), &
404	sk(k,j,i-1), u_comp, &
405	0.5, ddx )
406	ENDDO
407	degraded_l = .TRUE.
408
409	CASE ( 7 )
410
411	DO k = nzb_s_inner(j,i)+1, nzt
412	v_comp = v(k,j+1,i)-v_gtrans
413	flux_n(k) = v_comp * ( &
414	7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
415	- ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
416	diss_n(k) = -ABS( v_comp ) * ( &
417	3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
418	- ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
419	ENDDO
420
421	CASE ( 8 )
422
423	DO k = nzb_s_inner(j,i)+1, nzt
424	v_comp = v(k,j+1,i) - v_gtrans
425	flux_n(k) = v_comp * ( &
426	7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
427	- ( sk(k,j+2,i) + sk(k,j-1,i) ) ) * adv_sca_3
428	diss_n(k) = -ABS( v_comp ) * ( &
429	3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
430	- ( sk(k,j+2,i) - sk(k,j-1,i) ) ) * adv_sca_3
431	!
432	!-- Compute southside fluxes for the south boundary of PE domain
433	v_comp = v(k,j,i) - v_gtrans
434	swap_flux_y_local(k,tn) = v_comp * &
435	( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
436	swap_diss_y_local(k,tn) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), &
437	sk(k,j,i), sk(k,j-1,i), &
438	sk(k,j-1,i), v_comp, &
439	0.5, ddy )
440	ENDDO
441	degraded_s = .TRUE.
442
443	CASE DEFAULT
444
445	END SELECT
446
447	!
448	!-- Compute the crosswise 5th order fluxes at the outflow
449	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. &
450	boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN
451
452	DO k = nzb_s_inner(j,i)+1, nzt
453	v_comp = v(k,j+1,i) - v_gtrans
454	flux_n(k) = v_comp * ( &
455	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
456	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
457	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
458	diss_n(k) = -ABS( v_comp ) * ( &
459	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
460	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
461	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
462	ENDDO
463
464	ELSE
465
466	DO k = nzb_s_inner(j,i)+1, nzt
467	u_comp = u(k,j,i+1) - u_gtrans
468	flux_r(k) = u_comp * ( &
469	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
470	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
471	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
472	diss_r(k) = -ABS( u_comp ) * ( &
473	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
474	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
475	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
476	ENDDO
477
478	ENDIF
479
480	ELSE
481
482	!
483	!-- Compute the fifth order fluxes for the interior of PE domain.
484	DO k = nzb_u_inner(j,i)+1, nzt
485	u_comp = u(k,j,i+1) - u_gtrans
486	flux_r(k) = u_comp * ( &
487	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
488	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
489	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) * adv_sca_5
490	diss_r(k) = -ABS( u_comp ) * ( &
491	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
492	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
493	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) * adv_sca_5
494
495	v_comp = v(k,j+1,i) - v_gtrans
496	flux_n(k) = v_comp * ( &
497	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
498	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
499	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) * adv_sca_5
500	diss_n(k) = -ABS( v_comp ) * ( &
501	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
502	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
503	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5
504	ENDDO
505
506	ENDIF
507	!
508	!-- Compute left- and southside fluxes of the respective PE bounds.
509	IF ( j == nys .AND. .NOT. degraded_s ) THEN
510
511	DO k = nzb_s_inner(j,i)+1, nzt
512	v_comp = v(k,j,i) - v_gtrans
513	swap_flux_y_local(k,tn) = v_comp * ( &
514	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
515	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
516	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) &
517	) * adv_sca_5
518	swap_diss_y_local(k,tn) = -ABS( v_comp ) * ( &
519	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
520	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
521	+ sk(k,j+2,i) - sk(k,j-3,i) &
522	) * adv_sca_5
523	ENDDO
524
525	ENDIF
526
527	IF ( i == i_omp .AND. .NOT. degraded_l ) THEN
528
529	DO k = nzb_s_inner(j,i)+1, nzt
530	u_comp = u(k,j,i) - u_gtrans
531	swap_flux_x_local(k,j,tn) = u_comp * ( &
532	37.0 * ( sk(k,j,i) + sk(k,j,i-1) ) &
533	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
534	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) &
535	) * adv_sca_5
536	swap_diss_x_local(k,j,tn) = -ABS( u_comp ) * ( &
537	10.0 * ( sk(k,j,i) - sk(k,j,i-1) ) &
538	- 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
539	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) &
540	) * adv_sca_5
541	ENDDO
542
543	ENDIF
544
545	!
546	!-- Now compute the tendency terms for the horizontal parts
547	DO k = nzb_s_inner(j,i)+1, nzt
548
549	tend(k,j,i) = tend(k,j,i) - ( &
550	( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &
551	swap_diss_x_local(k,j,tn) ) * ddx &
552	+ ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &
553	swap_diss_y_local(k,tn) ) * ddy &
554	)
555
556	swap_flux_y_local(k,tn) = flux_n(k)
557	swap_diss_y_local(k,tn) = diss_n(k)
558	swap_flux_x_local(k,j,tn) = flux_r(k)
559	swap_diss_x_local(k,j,tn) = diss_r(k)
560
561	ENDDO
562
563	!
564	!-- Vertical advection, degradation of order near bottom and top.
565	!-- The fluxes flux_d and diss_d at the surface are 0. Due to later
566	!-- calculation of statistics the top flux at the surface should be 0.
567	flux_t(nzb_s_inner(j,i)) = 0.0
568	diss_t(nzb_s_inner(j,i)) = 0.0
569
570	!
571	!-- 2nd-order scheme (bottom)
572	k = nzb_s_inner(j,i)+1
573	flux_d = flux_t(k-1)
574	diss_d = diss_t(k-1)
575	flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
576
577	!
578	!-- sk(k,j,i) is referenced three times to avoid an access below surface
579	diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), sk(k,j,i), &
580	sk(k,j,i), w(k,j,i), 0.5, ddzw(k) )
581
582	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
583	* ddzw(k)
584	!
585	!-- WS3 as an intermediate step (bottom)
586	k = nzb_s_inner(j,i) + 2
587	flux_d = flux_t(k-1)
588	diss_d = diss_t(k-1)
589	flux_t(k) = w(k,j,i) * ( &
590	7.0 * ( sk(k+1,j,i) + sk(k,j,i) ) &
591	- ( sk(k+2,j,i) + sk(k-1,j,i) ) &
592	) * adv_sca_3
593	diss_t(k) = -ABS( w(k,j,i) ) * ( &
594	3.0 * ( sk(k+1,j,i) - sk(k,j,i) ) &
595	- ( sk(k+2,j,i) - sk(k-1,j,i) ) &
596	) * adv_sca_3
597
598	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
599	* ddzw(k)
600	!
601	!-- WS5
602	DO k = nzb_s_inner(j,i)+3, nzt-2
603
604	flux_d = flux_t(k-1)
605	diss_d = diss_t(k-1)
606	flux_t(k) = w(k,j,i) * ( &
607	37.0 * ( sk(k+1,j,i) + sk(k,j,i) ) &
608	- 8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) ) &
609	+ ( sk(k+3,j,i) + sk(k-2,j,i) ) &
610	) * adv_sca_5
611	diss_t(k) = -ABS( w(k,j,i) ) * ( &
612	10.0 * ( sk(k+1,j,i) - sk(k,j,i) )&
613	- 5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) )&
614	+ ( sk(k+3,j,i) - sk(k-2,j,i) )&
615	) * adv_sca_5
616
617	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
618	- ( flux_d + diss_d ) ) * ddzw(k)
619
620	ENDDO
621
622	!
623	!-- WS3 as an intermediate step (top)
624	k = nzt - 1
625	flux_d = flux_t(k-1)
626	diss_d = diss_t(k-1)
627	flux_t(k) = w(k,j,i) * ( &
628	7.0 * ( sk(k+1,j,i) + sk(k,j,i) ) &
629	- ( sk(k+2,j,i) + sk(k-1,j,i) ) &
630	) * adv_sca_3
631	diss_t(k) = -ABS( w(k,j,i) ) * ( &
632	3.0 * ( sk(k+1,j,i) - sk(k,j,i) ) &
633	- ( sk(k+2,j,i) - sk(k-1,j,i) ) &
634	) * adv_sca_3
635
636	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
637	* ddzw(k)
638	!
639	!-- 2nd-order scheme (top)
640	k = nzt
641	flux_d = flux_t(k-1)
642	diss_d = diss_t(k-1)
643	flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
644
645	!
646	!-- sk(k+1) is referenced two times to avoid a segmentation fault at top
647	diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), sk(k-1,j,i), &
648	sk(k-2,j,i), w(k,j,i), 0.5, ddzw(k) )
649
650	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) - flux_d - diss_d ) &
651	* ddzw(k)
652	!
653	!-- Evaluation of statistics
654	SELECT CASE ( sk_char )
655
656	CASE ( 'pt' )
657
658	DO k = nzb_s_inner(j,i), nzt
659	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) + &
660	( flux_t(k) + diss_t(k) ) &
661	* weight_substep(intermediate_timestep_count)
662	ENDDO
663
664	CASE ( 'sa' )
665
666	DO k = nzb_s_inner(j,i), nzt
667	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) + &
668	( flux_t(k) + diss_t(k) ) &
669	* weight_substep(intermediate_timestep_count)
670	ENDDO
671
672	CASE ( 'q' )
673
674	DO k = nzb_s_inner(j,i), nzt
675	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) + &
676	( flux_t(k) + diss_t(k) ) &
677	* weight_substep(intermediate_timestep_count)
678	ENDDO
679
680	END SELECT
681
682	END SUBROUTINE advec_s_ws_ij
683
684
685
686
687	!------------------------------------------------------------------------------!
688	! Advection of u-component - Call for grid point i,j
689	!------------------------------------------------------------------------------!
690	SUBROUTINE advec_u_ws_ij( i, j, i_omp, tn )
691
692	USE arrays_3d
693	USE constants
694	USE control_parameters
695	USE grid_variables
696	USE indices
697	USE statistics
698
699	IMPLICIT NONE
700
701	INTEGER :: i, i_omp, j, k, tn
702	LOGICAL :: degraded_l, degraded_s
703	REAL :: gu, gv, flux_d, diss_d, u_comp_l, v_comp, w_comp
704	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, &
705	flux_n, diss_n, u_comp
706
707	degraded_l = .FALSE.
708	degraded_s = .FALSE.
709
710	gu = 2.0 * u_gtrans
711	gv = 2.0 * v_gtrans
712
713	IF ( boundary_flags(j,i) /= 0 ) THEN
714	!
715	!-- Degrade the order for Dirichlet bc. at the outflow boundary
716	SELECT CASE ( boundary_flags(j,i) )
717
718	CASE ( 1 )
719	DO k = nzb_u_inner(j,i)+1, nzt
720	u_comp(k) = u(k,j,i+1) + u(k,j,i)
721	flux_r(k) = ( u_comp(k) - gu ) * ( &
722	7.0 * ( u(k,j,i+1) + u(k,j,i) ) &
723	- ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
724	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
725	3.0 * ( u(k,j,i+1) - u(k,j,i) ) &
726	- ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
727	ENDDO
728
729	CASE ( 2 )
730	DO k = nzb_u_inner(j,i)+1, nzt
731	u_comp(k) = u(k,j,i+1) + u(k,j,i)
732	flux_r(k) = ( u_comp(k) - gu ) * ( &
733	u(k,j,i+1) + u(k,j,i) ) * 0.25
734	diss_r(k) = diss_2nd( u(k,j,i+1) ,u(k,j,i+1), u(k,j,i), &
735	u(k,j,i-1), u(k,j,i-2), u_comp(k), &
736	0.25, ddx )
737	ENDDO
738
739	CASE ( 3 )
740	DO k = nzb_u_inner(j,i)+1, nzt
741	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
742	flux_n(k) = v_comp * ( &
743	7.0 * ( u(k,j+1,i) + u(k,j,i) ) &
744	- ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
745	diss_n(k) = - ABS( v_comp ) * ( &
746	3.0 * ( u(k,j+1,i) - u(k,j,i) ) &
747	- ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
748	ENDDO
749
750	CASE ( 4 )
751	DO k = nzb_u_inner(j,i)+1, nzt
752	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
753	flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
754	diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), u(k,j,i), &
755	u(k,j-1,i), u(k,j-2,i), v_comp, &
756	0.25, ddy )
757	ENDDO
758
759	CASE ( 5 )
760	DO k = nzb_u_inner(j,i)+1, nzt
761	!
762	!-- Compute leftside fluxes for the left boundary of PE domain
763	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
764	flux_l_u(k,j,tn) = u_comp(k) * ( u(k,j,i) + u(k,j,i-1) ) * 0.25
765	diss_l_u(k,j,tn) = diss_2nd( u(k,j,i+2), u(k,j,i+1), u(k,j,i), &
766	u(k,j,i-1), u(k,j,i-1), u_comp(k),&
767	0.25, ddx )
768
769	u_comp(k) = u(k,j,i+1) + u(k,j,i)
770	flux_r(k) = ( u_comp(k) - gu ) * ( &
771	7.0 * ( u(k,j,i+1) + u(k,j,i) ) &
772	- ( u(k,j,i+2) + u(k,j,i-1) ) ) * adv_mom_3
773	diss_r(k) = - ABS( u_comp(k) -gu ) * ( &
774	3.0 * ( u(k,j,i+1) - u(k,j,i) ) &
775	- ( u(k,j,i+2) - u(k,j,i-1) ) ) * adv_mom_3
776	ENDDO
777	degraded_l = .TRUE.
778
779	CASE ( 7 )
780	DO k = nzb_u_inner(j,i)+1, nzt
781	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
782	flux_n(k) = v_comp * ( &
783	7.0 * ( u(k,j+1,i) + u(k,j,i) ) &
784	- ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
785	diss_n(k) = - ABS( v_comp ) * ( &
786	3.0 * ( u(k,j+1,i) - u(k,j,i) ) &
787	- ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
788	ENDDO
789
790	CASE ( 8 )
791	DO k = nzb_u_inner(j,i)+1, nzt
792	!
793	!-- Compute southside fluxes for the south boundary of PE domain
794	v_comp = v(k,j,i) + v(k,j,i-1) - gv
795	flux_s_u(k,tn) = v_comp * ( u(k,j,i) + u(k,j-1,i) ) * 0.25
796	diss_s_u(k,tn) = diss_2nd( u(k,j+2,i), u(k,j+1,i), u(k,j,i), &
797	u(k,j-1,i), u(k,j-1,i), v_comp, &
798	0.25, ddy )
799
800	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
801	flux_n(k) = v_comp * ( &
802	7.0 * ( u(k,j+1,i) + u(k,j,i) ) &
803	- ( u(k,j+2,i) + u(k,j-1,i) ) ) * adv_mom_3
804	diss_n(k) = - ABS( v_comp ) * ( &
805	3.0 * ( u(k,j+1,i) - u(k,j,i) ) &
806	- ( u(k,j+2,i) - u(k,j-1,i) ) ) * adv_mom_3
807	ENDDO
808	degraded_s = .TRUE.
809
810	CASE DEFAULT
811
812	END SELECT
813	!
814	!-- Compute the crosswise 5th order fluxes at the outflow
815	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. &
816	boundary_flags(j,i) == 5 ) THEN
817
818	DO k = nzb_u_inner(j,i)+1, nzt
819	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
820	flux_n(k) = v_comp * ( &
821	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
822	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
823	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
824	diss_n(k) = - ABS ( v_comp ) * ( &
825	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
826	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
827	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
828	ENDDO
829
830	ELSE
831
832	DO k = nzb_u_inner(j,i)+1, nzt
833	u_comp(k) = u(k,j,i+1) + u(k,j,i)
834	flux_r(k) = ( u_comp(k) - gu ) * ( &
835	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
836	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
837	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
838	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
839	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
840	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
841	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
842	ENDDO
843
844	ENDIF
845
846	ELSE
847	!
848	!-- Compute the fifth order fluxes for the interior of PE domain.
849	DO k = nzb_u_inner(j,i)+1, nzt
850	u_comp(k) = u(k,j,i+1) + u(k,j,i)
851	flux_r(k) = ( u_comp(k) - gu ) * ( &
852	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
853	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
854	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5
855	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
856	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
857	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
858	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
859
860	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
861	flux_n(k) = v_comp * ( &
862	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
863	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
864	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
865	diss_n(k) = - ABS( v_comp ) * ( &
866	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
867	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
868	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
869	ENDDO
870
871	ENDIF
872	!
873	!-- Compute left- and southside fluxes for the respective boundary of PE
874	IF ( j == nys .AND. ( .NOT. degraded_s ) ) THEN
875
876	DO k = nzb_u_inner(j,i)+1, nzt
877	v_comp = v(k,j,i) + v(k,j,i-1) - gv
878	flux_s_u(k,tn) = v_comp * ( &
879	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
880	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
881	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) * adv_mom_5
882	diss_s_u(k,tn) = - ABS(v_comp) * ( &
883	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
884	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
885	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) * adv_mom_5
886	ENDDO
887
888	ENDIF
889
890	IF ( i == i_omp .AND. ( .NOT. degraded_l ) ) THEN
891
892	DO k = nzb_u_inner(j,i)+1, nzt
893	u_comp_l = u(k,j,i)+u(k,j,i-1)-gu
894	flux_l_u(k,j,tn) = u_comp_l * ( &
895	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
896	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
897	+ ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5
898	diss_l_u(k,j,tn) = - ABS(u_comp_l) * ( &
899	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
900	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
901	+ ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5
902	ENDDO
903
904	ENDIF
905	!
906	!-- Now compute the tendency terms for the horizontal parts.
907	DO k = nzb_u_inner(j,i)+1, nzt
908	tend(k,j,i) = tend(k,j,i) - ( &
909	( flux_r(k) + diss_r(k) &
910	- flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &
911	+ ( flux_n(k) + diss_n(k) &
912	- flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy )
913
914	flux_l_u(k,j,tn) = flux_r(k)
915	diss_l_u(k,j,tn) = diss_r(k)
916	flux_s_u(k,tn) = flux_n(k)
917	diss_s_u(k,tn) = diss_n(k)
918	!
919	!-- Statistical Evaluation of u'u'. The factor has to be applied for
920	!-- right evaluation when gallilei_trans = .T. .
921	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
922	+ ( flux_r(k) * &
923	( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
924	/ ( u_comp(k) - gu + 1.0E-20 ) &
925	+ diss_r(k) * &
926	ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
927	/ ( ABS( u_comp(k) - gu ) + 1.0E-20 ) ) &
928	* weight_substep(intermediate_timestep_count)
929	ENDDO
930	sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
931
932
933	!
934	!-- Vertical advection, degradation of order near surface and top.
935	!-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
936	!-- statistical evaluation the top flux at the surface should be 0
937	flux_t(nzb_u_inner(j,i)) = 0. !statistical reasons
938	diss_t(nzb_u_inner(j,i)) = 0.
939	!
940	!-- 2nd order scheme (bottom)
941	k = nzb_u_inner(j,i)+1
942	flux_d = flux_t(k-1)
943	diss_d = diss_t(k-1)
944	w_comp = w(k,j,i) + w(k,j,i-1)
945	flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) *0.25
946	diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), 0., 0., &
947	w_comp, 0.25, ddzw(k) )
948
949	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
950	- flux_d - diss_d ) * ddzw(k)
951	!
952	!-- WS3 as an intermediate step (bottom)
953	k = nzb_u_inner(j,i)+2
954	flux_d = flux_t(k-1)
955	diss_d = diss_t(k-1)
956	w_comp = w(k,j,i) + w(k,j,i-1)
957	flux_t(k) = w_comp * ( &
958	7.0 * ( u(k+1,j,i) + u(k,j,i) ) &
959	- ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
960	diss_t(k) = -ABS( w_comp) * ( &
961	3.0 * ( u(k+1,j,i) - u(k,j,i) ) &
962	- ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
963
964	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
965	- flux_d - diss_d ) * ddzw(k)
966	!
967	!-- WS5
968	DO k = nzb_u_inner(j,i)+3, nzt-2
969	flux_d = flux_t(k-1)
970	diss_d = diss_t(k-1)
971	w_comp = w(k,j,i) + w(k,j,i-1)
972	flux_t(k) = w_comp * ( &
973	37.0 * ( u(k+1,j,i) + u(k,j,i) ) &
974	- 8.0 * ( u(k+2,j,i) + u(k-1,j,i) ) &
975	+ ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
976	diss_t(k) = - ABS( w_comp ) * ( &
977	10.0 * ( u(k+1,j,i) - u(k,j,i) ) &
978	- 5.0 * ( u(k+2,j,i) - u(k-1,j,i) ) &
979	+ ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
980
981	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
982	- flux_d - diss_d ) * ddzw(k)
983	ENDDO
984	!
985	!-- WS3 as an intermediate step (top)
986	k = nzt - 1
987	flux_d = flux_t(k-1)
988	diss_d = diss_t(k-1)
989	w_comp = w(k,j,i) + w(k,j,i-1)
990	flux_t(k) = w_comp * ( &
991	7.0 * ( u(k+1,j,i) + u(k,j,i) ) &
992	- ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
993	diss_t(k) = - ABS( w_comp ) * ( &
994	3.0 * ( u(k+1,j,i) - u(k,j,i) ) &
995	- ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
996
997	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
998	- flux_d - diss_d ) * ddzw(k)
999
1000	!
1001	!-- 2nd order scheme (top)
1002	k = nzt
1003	flux_d = flux_t(k-1)
1004	diss_d = diss_t(k-1)
1005	w_comp = w(k,j,i)+w(k,j,i-1)
1006	flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
1007	diss_t(k) = diss_2nd( u(k+1,j,i), u(k+1,j,i), u(k,j,i), u(k-1,j,i), &
1008	u(k-2,j,i), w_comp, 0.25, ddzw(k) )
1009
1010	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1011	- flux_d - diss_d ) * ddzw(k)
1012
1013	!
1014	!-- sum up the vertical momentum fluxes
1015	DO k = nzb_u_inner(j,i), nzt
1016	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
1017	+ ( flux_t(k) + diss_t(k) ) &
1018	* weight_substep(intermediate_timestep_count)
1019	ENDDO
1020
1021
1022	END SUBROUTINE advec_u_ws_ij
1023
1024
1025
1026
1027	!-----------------------------------------------------------------------------!
1028	! Advection of v-component - Call for grid point i,j
1029	!-----------------------------------------------------------------------------!
1030	SUBROUTINE advec_v_ws_ij( i, j, i_omp, tn )
1031
1032	USE arrays_3d
1033	USE constants
1034	USE control_parameters
1035	USE grid_variables
1036	USE indices
1037	USE statistics
1038
1039	IMPLICIT NONE
1040
1041	INTEGER :: i, i_omp, j, k, tn
1042	LOGICAL :: degraded_l, degraded_s
1043	REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp_l, w_comp
1044	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, &
1045	diss_n, flux_r, diss_r, v_comp
1046
1047	degraded_l = .FALSE.
1048	degraded_s = .FALSE.
1049
1050	gu = 2.0 * u_gtrans
1051	gv = 2.0 * v_gtrans
1052
1053	IF ( boundary_flags(j,i) /= 0 ) THEN
1054	!
1055	!-- Degrade the order for Dirichlet bc. at the outflow boundary
1056	SELECT CASE ( boundary_flags(j,i) )
1057
1058	CASE ( 1 )
1059	DO k = nzb_v_inner(j,i)+1, nzt
1060	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1061	flux_r(k) = u_comp * ( &
1062	7.0 * (v(k,j,i+1) + v(k,j,i) ) &
1063	- ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1064	diss_r(k) = - ABS( u_comp ) * ( &
1065	3.0 * ( v(k,j,i+1) - v(k,j,i) ) &
1066	- ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
1067	ENDDO
1068
1069	CASE ( 2 )
1070	DO k = nzb_v_inner(j,i)+1, nzt
1071	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1072	flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25
1073	diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), v(k,j,i), &
1074	v(k,j,i-1), v(k,j,i-2), u_comp, &
1075	0.25, ddx )
1076	ENDDO
1077
1078	CASE ( 3 )
1079	DO k = nzb_v_inner(j,i)+1, nzt
1080	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1081	flux_n(k) = ( v_comp(k)- gv ) * ( &
1082	7.0 * ( v(k,j+1,i) + v(k,j,i) ) &
1083	- ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
1084	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
1085	3.0 * ( v(k,j+1,i) - v(k,j,i) ) &
1086	- ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
1087	ENDDO
1088
1089	CASE ( 4 )
1090	DO k = nzb_v_inner(j,i)+1, nzt
1091	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1092	flux_n(k) = ( v_comp(k) - gv ) * &
1093	( v(k,j+1,i) + v(k,j,i) ) * 0.25
1094	diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), v(k,j,i), &
1095	v(k,j-1,i), v(k,j-2,i), v_comp(k), &
1096	0.25, ddy )
1097	ENDDO
1098
1099	CASE ( 5 )
1100	DO k = nzb_w_inner(j,i)+1, nzt
1101	u_comp = u(k,j-1,i) + u(k,j,i) - gu
1102	flux_r(k) = u_comp * ( &
1103	7.0 * ( v(k,j,i+1) + v(k,j,i) ) &
1104	- ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1105	diss_r(k) = - ABS( u_comp ) * ( &
1106	3.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1107	- ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
1108	ENDDO
1109
1110	CASE ( 6 )
1111	DO k = nzb_v_inner(j,i)+1, nzt
1112	u_comp = u(k,j-1,i) + u(k,j,i) - gu
1113	flux_l_v(k,j,tn) = u_comp * ( v(k,j,i) + v(k,j,i-1) ) * 0.25
1114	diss_l_v(k,j,tn) = diss_2nd( v(k,j,i+2), v(k,j,i+1), v(k,j,i),&
1115	v(k,j,i-1), v(k,j,i-1), u_comp, &
1116	0.25, ddx )
1117
1118	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1119	flux_r(k) = u_comp * ( &
1120	7.0 * ( v(k,j,i+1) + v(k,j,i) ) &
1121	- ( v(k,j,i+2) + v(k,j,i-1) ) ) * adv_mom_3
1122	diss_r(k) = - ABS( u_comp ) * ( &
1123	3.0 * ( v(k,j,i+1) - v(k,j,i) ) &
1124	- ( v(k,j,i+2) - v(k,j,i-1) ) ) * adv_mom_3
1125	ENDDO
1126	degraded_l = .TRUE.
1127
1128	CASE ( 7 )
1129	DO k = nzb_v_inner(j,i)+1, nzt
1130	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
1131	flux_s_v(k,tn) = v_comp(k) * ( v(k,j,i) + v(k,j-1,i) ) * 0.25
1132	diss_s_v(k,tn) = diss_2nd( v(k,j+2,i), v(k,j+1,i), v(k,j,i), &
1133	v(k,j-1,i), v(k,j-1,i), v_comp(k), &
1134	0.25, ddy )
1135
1136	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1137	flux_n(k) = ( v_comp(k) - gv ) * ( &
1138	7.0 * ( v(k,j+1,i) + v(k,j,i) ) &
1139	- ( v(k,j+2,i) + v(k,j-1,i) ) ) * adv_mom_3
1140	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
1141	3.0 * ( v(k,j+1,i) - v(k,j,i) ) &
1142	- ( v(k,j+2,i) - v(k,j-1,i) ) ) * adv_mom_3
1143	ENDDO
1144	degraded_s = .TRUE.
1145
1146	CASE DEFAULT
1147
1148	END SELECT
1149	!
1150	!-- Compute the crosswise 5th order fluxes at the outflow
1151	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. &
1152	boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN
1153
1154	DO k = nzb_v_inner(j,i)+1, nzt
1155	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1156	flux_n(k) = ( v_comp(k) - gv ) * ( &
1157	37.0 * ( v(k,j+1,i) + v(k,j,i) ) &
1158	- 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) &
1159	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1160	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
1161	10.0 * ( v(k,j+1,i) - v(k,j,i) ) &
1162	- 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) &
1163	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1164	ENDDO
1165
1166	ELSE
1167
1168	DO k = nzb_v_inner(j,i)+1, nzt
1169	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1170	flux_r(k) = u_comp * ( &
1171	37.0 * ( v(k,j,i+1) + v(k,j,i) ) &
1172	- 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) &
1173	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1174	diss_r(k) = - ABS( u_comp ) * ( &
1175	10.0 * ( v(k,j,i+1) - v(k,j,i) ) &
1176	- 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) &
1177	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1178	ENDDO
1179
1180	ENDIF
1181
1182	ELSE
1183	!
1184	!-- Compute the fifth order fluxes for the interior of PE domain.
1185	DO k = nzb_v_inner(j,i)+1, nzt
1186	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
1187	flux_r(k) = u_comp * ( &
1188	37.0 * ( v(k,j,i+1) + v(k,j,i) ) &
1189	- 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) &
1190	+ ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5
1191	diss_r(k) = - ABS( u_comp ) * ( &
1192	10.0 * ( v(k,j,i+1) - v(k,j,i) ) &
1193	- 5.0 * ( v(k,j,i+2) - v(k,j,i-1) ) &
1194	+ ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5
1195
1196	v_comp(k) = v(k,j+1,i) + v(k,j,i)
1197	flux_n(k) = ( v_comp(k) - gv ) * ( &
1198	37.0 * ( v(k,j+1,i) + v(k,j,i) ) &
1199	- 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) &
1200	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5
1201	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
1202	10.0 * ( v(k,j+1,i) - v(k,j,i) ) &
1203	- 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) &
1204	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5
1205	ENDDO
1206
1207	ENDIF
1208	!
1209	!-- Compute left- and southside fluxes for the respective boundary
1210	IF ( i == i_omp .AND. .NOT. degraded_l ) THEN
1211
1212	DO k = nzb_v_inner(j,i)+1, nzt
1213	u_comp = u(k,j-1,i) + u(k,j,i) - gu
1214	flux_l_v(k,j,tn) = u_comp * ( &
1215	37.0 * ( v(k,j,i) + v(k,j,i-1) ) &
1216	- 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) &
1217	+ ( v(k,j,i+2) + v(k,j,i-3) ) ) * adv_mom_5
1218	diss_l_v(k,j,tn) = - ABS( u_comp ) * ( &
1219	10.0 * ( v(k,j,i) - v(k,j,i-1) ) &
1220	- 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) &
1221	+ ( v(k,j,i+2) - v(k,j,i-3) ) ) * adv_mom_5
1222	ENDDO
1223
1224	ENDIF
1225
1226	IF ( j == nysv .AND. .NOT. degraded_s ) THEN
1227
1228	DO k = nzb_v_inner(j,i)+1, nzt
1229	v_comp_l = v(k,j,i) + v(k,j-1,i) - gv
1230	flux_s_v(k,tn) = v_comp_l * ( &
1231	37.0 * ( v(k,j,i) + v(k,j-1,i) ) &
1232	- 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) &
1233	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) * adv_mom_5
1234	diss_s_v(k,tn) = - ABS( v_comp_l ) * ( &
1235	10.0 * ( v(k,j,i) - v(k,j-1,i) ) &
1236	- 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) &
1237	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) * adv_mom_5
1238	ENDDO
1239
1240	ENDIF
1241	!
1242	!-- Now compute the tendency terms for the horizontal parts.
1243	DO k = nzb_v_inner(j,i)+1, nzt
1244	tend(k,j,i) = tend(k,j,i) - ( &
1245	( flux_r(k) + diss_r(k) &
1246	- flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &
1247	+ ( flux_n(k) + diss_n(k) &
1248	- flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy )
1249
1250	flux_l_v(k,j,tn) = flux_r(k)
1251	diss_l_v(k,j,tn) = diss_r(k)
1252	flux_s_v(k,tn) = flux_n(k)
1253	diss_s_v(k,tn) = diss_n(k)
1254
1255	!
1256	!-- Statistical Evaluation of v'v'. The factor has to be applied for
1257	!-- right evaluation when gallilei_trans = .T. .
1258
1259	sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &
1260	+ ( flux_n(k) &
1261	* ( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
1262	/ ( v_comp(k) - gv + 1.0E-20 ) &
1263	+ diss_n(k) &
1264	* ABS( v_comp(k) - 2.0 * hom(k,1,2,0) ) &
1265	/ ( ABS( v_comp(k) - gv ) +1.0E-20 ) ) &
1266	* weight_substep(intermediate_timestep_count)
1267
1268	ENDDO
1269	sums_vs2_ws_l(nzb_v_inner(j,i),tn) = sums_vs2_ws_l(nzb_v_inner(j,i)+1,tn)
1270
1271	!
1272	!-- Vertical advection, degradation of order near surface and top.
1273	!-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1274	!-- statistical evaluation the top flux at the surface should be 0
1275	flux_t(nzb_v_inner(j,i)) = 0.0 !statistical reasons
1276	diss_t(nzb_v_inner(j,i)) = 0.0
1277	!
1278	!-- 2nd order scheme (bottom)
1279	k = nzb_v_inner(j,i)+1
1280	flux_d = flux_t(k-1)
1281	diss_d = diss_t(k-1)
1282	w_comp = w(k,j-1,i) + w(k,j,i)
1283	flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
1284	diss_t(k) = diss_2nd( v(k+2,j,i), v(k+1,j,i), v(k,j,i), 0.0, 0.0, &
1285	w_comp, 0.25, ddzw(k) )
1286
1287	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1288	- flux_d - diss_d ) * ddzw(k)
1289
1290	!
1291	!-- WS3 as an intermediate step (bottom)
1292	k = nzb_v_inner(j,i)+2
1293	flux_d = flux_t(k-1)
1294	diss_d = diss_t(k-1)
1295	w_comp = w(k,j-1,i) + w(k,j,i)
1296	flux_t(k) = w_comp * ( &
1297	7.0 * ( v(k+1,j,i) + v(k,j,i) ) &
1298	- ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
1299	diss_t(k) = - ABS( w_comp ) * ( &
1300	3.0 * ( v(k+1,j,i) - v(k,j,i) ) &
1301	- ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
1302
1303	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1304	- flux_d - diss_d ) * ddzw(k)
1305	!
1306	!-- WS5
1307	DO k = nzb_v_inner(j,i)+3, nzt-2
1308	flux_d = flux_t(k-1)
1309	diss_d = diss_t(k-1)
1310	w_comp = w(k,j-1,i) + w(k,j,i)
1311	flux_t(k) = w_comp * ( &
1312	37.0 * ( v(k+1,j,i) + v(k,j,i) ) &
1313	- 8.0 * ( v(k+2,j,i) + v(k-1,j,i) ) &
1314	+ ( v(k+3,j,i) + v(k-2,j,i) ) ) * adv_mom_5
1315	diss_t(k) = - ABS( w_comp ) * ( &
1316	10.0 * ( v(k+1,j,i) - v(k,j,i) ) &
1317	- 5.0 * ( v(k+2,j,i) - v(k-1,j,i) ) &
1318	+ ( v(k+3,j,i) - v(k-2,j,i) ) ) * adv_mom_5
1319
1320	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1321	- flux_d - diss_d ) * ddzw(k)
1322	ENDDO
1323	!
1324	!-- WS3 as an intermediate step (top)
1325	k = nzt - 1
1326	flux_d = flux_t(k-1)
1327	diss_d = diss_t(k-1)
1328	w_comp = w(k,j-1,i) + w(k,j,i)
1329	flux_t(k) = w_comp * ( &
1330	7.0 * ( v(k+1,j,i) + v(k,j,i) ) &
1331	- ( v(k+2,j,i) + v(k-1,j,i) ) ) * adv_mom_3
1332	diss_t(k) = - ABS( w_comp ) * ( &
1333	3.0 * ( v(k+1,j,i) - v(k,j,i) ) &
1334	- ( v(k+2,j,i) - v(k-1,j,i) ) ) * adv_mom_3
1335	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1336	- flux_d - diss_d ) * ddzw(k)
1337	!
1338	!-- 2nd order scheme (top)
1339	k = nzt
1340	flux_d = flux_t(k-1)
1341	diss_d = diss_t(k-1)
1342	w_comp = w(k,j-1,i)+w(k,j,i)
1343	flux_t(k) = w_comp * ( v(k+1,j,i) + v(k,j,i) ) * 0.25
1344	diss_t(k) = diss_2nd( v(k+1,j,i), v(k+1,j,i), v(k,j,i), v(k-1,j,i), &
1345	v(k-2,j,i), w_comp, 0.25, ddzw(k) )
1346
1347	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1348	- flux_d - diss_d ) * ddzw(k)
1349
1350	DO k = nzb_v_inner(j,i), nzt
1351	sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &
1352	+ ( flux_t(k) + diss_t(k) ) &
1353	* weight_substep(intermediate_timestep_count)
1354	ENDDO
1355
1356	END SUBROUTINE advec_v_ws_ij
1357
1358
1359
1360	!------------------------------------------------------------------------------!
1361	! Advection of w-component - Call for grid point i,j
1362	!------------------------------------------------------------------------------!
1363	SUBROUTINE advec_w_ws_ij( i, j, i_omp, tn )
1364
1365	USE arrays_3d
1366	USE constants
1367	USE control_parameters
1368	USE grid_variables
1369	USE indices
1370	USE statistics
1371
1372	IMPLICIT NONE
1373
1374	INTEGER :: i, i_omp, j, k, tn
1375	LOGICAL :: degraded_l, degraded_s
1376	REAL :: gu, gv, flux_d, diss_d, u_comp, v_comp, w_comp
1377	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n, &
1378	diss_n
1379
1380	degraded_l = .FALSE.
1381	degraded_s = .FALSE.
1382
1383	gu = 2.0 * u_gtrans
1384	gv = 2.0 * v_gtrans
1385
1386
1387	IF ( boundary_flags(j,i) /= 0 ) THEN
1388	!
1389	!-- Degrade the order for Dirichlet bc. at the outflow boundary
1390	SELECT CASE ( boundary_flags(j,i) )
1391
1392	CASE ( 1 )
1393	DO k = nzb_w_inner(j,i)+1, nzt
1394	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1395	flux_r(k) = u_comp * ( &
1396	7.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1397	- ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1398	diss_r(k) = -ABS( u_comp ) * ( &
1399	3.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1400	- ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1401	ENDDO
1402
1403	CASE ( 2 )
1404	DO k = nzb_w_inner(j,i)+1, nzt
1405	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1406	flux_r(k) = u_comp * ( w(k,j,i+1) + w(k,j,i) ) * 0.25
1407	diss_r(k) = diss_2nd( w(k,j,i+1), w(k,j,i+1), w(k,j,i), &
1408	w(k,j,i-1), w(k,j,i-2), u_comp, &
1409	0.25, ddx )
1410	ENDDO
1411
1412	CASE ( 3 )
1413	DO k = nzb_w_inner(j,i)+1, nzt
1414	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1415	flux_n(k) = v_comp * ( &
1416	7.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1417	- ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1418	diss_n(k) = -ABS( v_comp ) * ( &
1419	3.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1420	- ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1421	ENDDO
1422
1423	CASE ( 4 )
1424	DO k = nzb_w_inner(j,i)+1, nzt
1425	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1426	flux_n(k) = v_comp * ( w(k,j+1,i) + w(k,j,i) ) * 0.25
1427	diss_n(k) = diss_2nd( w(k,j+1,i), w(k,j+1,i), w(k,j,i), &
1428	w(k,j-1,i), w(k,j-2,i), v_comp, &
1429	0.25, ddy )
1430	ENDDO
1431
1432	CASE ( 5 )
1433	DO k = nzb_w_inner(j,i)+1, nzt
1434	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1435	flux_r(k) = u_comp * ( &
1436	7.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1437	- ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1438	diss_r(k) = - ABS( u_comp ) * ( &
1439	3.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1440	- ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1441	ENDDO
1442
1443	CASE ( 6 )
1444	DO k = nzb_w_inner(j,i)+1, nzt
1445	!
1446	!-- Compute leftside fluxes for the left boundary of PE domain
1447	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1448	flux_r(k) = u_comp *( &
1449	7.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1450	- ( w(k,j,i+2) + w(k,j,i-1) ) ) * adv_mom_3
1451	diss_r(k) = - ABS( u_comp ) * ( &
1452	3.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1453	- ( w(k,j,i+2) - w(k,j,i-1) ) ) * adv_mom_3
1454
1455	u_comp = u(k+1,j,i) + u(k,j,i) - gu
1456	flux_l_w(k,j,tn) = u_comp * ( w(k,j,i) + w(k,j,i-1) ) * 0.25
1457	diss_l_w(k,j,tn) = diss_2nd( w(k,j,i+2), w(k,j,i+1), w(k,j,i), &
1458	w(k,j,i-1), w(k,j,i-1), u_comp, &
1459	0.25, ddx )
1460	ENDDO
1461	degraded_l = .TRUE.
1462
1463	CASE ( 7 )
1464	DO k = nzb_w_inner(j,i)+1, nzt
1465	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1466	flux_n(k) = v_comp *( &
1467	7.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1468	- ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1469	diss_n(k) = - ABS( v_comp ) * ( &
1470	3.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1471	- ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1472	ENDDO
1473
1474	CASE ( 8 )
1475	DO k = nzb_w_inner(j,i)+1, nzt
1476	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1477	flux_n(k) = v_comp * ( &
1478	7.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1479	- ( w(k,j+2,i) + w(k,j-1,i) ) ) * adv_mom_3
1480	diss_n(k) = - ABS( v_comp ) * ( &
1481	3.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1482	- ( w(k,j+2,i) - w(k,j-1,i) ) ) * adv_mom_3
1483	!
1484	!-- Compute southside fluxes for the south boundary of PE domain
1485	v_comp = v(k+1,j,i) + v(k,j,i) - gv
1486	flux_s_w(k,tn) = v_comp * ( w(k,j,i) + w(k,j-1,i) ) * 0.25
1487	diss_s_w(k,tn) = diss_2nd( w(k,j+2,i), w(k,j+1,i), w(k,j,i), &
1488	w(k,j-1,i), w(k,j-1,i), v_comp, &
1489	0.25, ddy )
1490	ENDDO
1491	degraded_s = .TRUE.
1492
1493	CASE DEFAULT
1494
1495	END SELECT
1496	!
1497	!-- Compute the crosswise 5th order fluxes at the outflow
1498	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR. &
1499	boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) THEN
1500
1501	DO k = nzb_w_inner(j,i)+1, nzt
1502	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1503	flux_n(k) = v_comp * ( &
1504	37.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1505	- 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) &
1506	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1507	diss_n(k) = - ABS( v_comp ) * ( &
1508	10.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1509	- 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) &
1510	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1511	ENDDO
1512
1513	ELSE
1514
1515	DO k = nzb_w_inner(j,i)+1, nzt
1516	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1517	flux_r(k) = u_comp * ( &
1518	37.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1519	- 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) &
1520	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1521	diss_r(k) = - ABS( u_comp ) * ( &
1522	10.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1523	- 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) &
1524	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1525	ENDDO
1526
1527	ENDIF
1528
1529	ELSE
1530	!
1531	!-- Compute the fifth order fluxes for the interior of PE domain.
1532	DO k = nzb_w_inner(j,i)+1, nzt
1533	u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu
1534	flux_r(k) = u_comp * ( &
1535	37.0 * ( w(k,j,i+1) + w(k,j,i) ) &
1536	- 8.0 * ( w(k,j,i+2) + w(k,j,i-1) ) &
1537	+ ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5
1538	diss_r(k) = - ABS( u_comp ) * ( &
1539	10.0 * ( w(k,j,i+1) - w(k,j,i) ) &
1540	- 5.0 * ( w(k,j,i+2) - w(k,j,i-1) ) &
1541	+ ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5
1542
1543	v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv
1544	flux_n(k) = v_comp * ( &
1545	37.0 * ( w(k,j+1,i) + w(k,j,i) ) &
1546	- 8.0 * ( w(k,j+2,i) + w(k,j-1,i) ) &
1547	+ ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5
1548	diss_n(k) = - ABS( v_comp ) * ( &
1549	10.0 * ( w(k,j+1,i) - w(k,j,i) ) &
1550	- 5.0 * ( w(k,j+2,i) - w(k,j-1,i) ) &
1551	+ ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5
1552	ENDDO
1553
1554	ENDIF
1555	!
1556	!-- Compute left- and southside fluxes for the respective boundary
1557	IF ( j == nys .AND. .NOT. degraded_s ) THEN
1558
1559	DO k = nzb_w_inner(j,i)+1, nzt
1560	v_comp = v(k+1,j,i) + v(k,j,i) - gv
1561	flux_s_w(k,tn) = v_comp * ( &
1562	37.0 * ( w(k,j,i) + w(k,j-1,i) ) &
1563	- 8.0 * ( w(k,j+1,i) +w(k,j-2,i) ) &
1564	+ ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5
1565	diss_s_w(k,tn) = - ABS( v_comp ) * ( &
1566	10.0 * ( w(k,j,i) - w(k,j-1,i) ) &
1567	- 5.0 * ( w(k,j+1,i) - w(k,j-2,i) ) &
1568	+ ( w(k,j+2,i) - w(k,j-3,i) ) ) * adv_mom_5
1569	ENDDO
1570
1571	ENDIF
1572
1573	IF ( i == i_omp .AND. .NOT. degraded_l ) THEN
1574
1575	DO k = nzb_w_inner(j,i)+1, nzt
1576	u_comp = u(k+1,j,i) + u(k,j,i) - gu
1577	flux_l_w(k,j,tn) = u_comp * ( &
1578	37.0 * ( w(k,j,i) + w(k,j,i-1) ) &
1579	- 8.0 * ( w(k,j,i+1) + w(k,j,i-2) ) &
1580	+ ( w(k,j,i+2) + w(k,j,i-3) ) ) * adv_mom_5
1581	diss_l_w(k,j,tn) = - ABS( u_comp ) * ( &
1582	10.0 * ( w(k,j,i) - w(k,j,i-1) ) &
1583	- 5.0 * ( w(k,j,i+1) - w(k,j,i-2) ) &
1584	+ ( w(k,j,i+2) - w(k,j,i-3) ) ) * adv_mom_5
1585	ENDDO
1586
1587	ENDIF
1588
1589	!
1590	!-- Now compute the tendency terms for the horizontal parts.
1591	DO k = nzb_w_inner(j,i)+1, nzt
1592	tend(k,j,i) = tend(k,j,i) - ( &
1593	( flux_r(k) + diss_r(k) &
1594	- flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &
1595	+ ( flux_n(k) + diss_n(k) &
1596	- flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy )
1597
1598	flux_l_w(k,j,tn) = flux_r(k)
1599	diss_l_w(k,j,tn) = diss_r(k)
1600	flux_s_w(k,tn) = flux_n(k)
1601	diss_s_w(k,tn) = diss_n(k)
1602	ENDDO
1603
1604	!
1605	!-- Vertical advection, degradation of order near surface and top.
1606	!-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1607	!-- statistical evaluation the top flux at the surface should be 0
1608	flux_t(nzb_w_inner(j,i)) = 0.0 !statistical reasons
1609	diss_t(nzb_w_inner(j,i)) = 0.0
1610	!
1611	!-- 2nd order scheme (bottom)
1612	k = nzb_w_inner(j,i)+1
1613	flux_d = flux_t(k-1)
1614	diss_d = diss_t(k-1)
1615	w_comp = w(k+1,j,i) + w(k,j,i)
1616	flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
1617	diss_t(k) = diss_2nd( w(k+2,j,i), w(k+1,j,i), w(k,j,i), 0.0, 0.0, &
1618	w_comp, 0.25, ddzu(k+1) )
1619
1620	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1621	- flux_d - diss_d ) * ddzu(k+1)
1622	!
1623	!-- WS3 as an intermediate step (bottom)
1624	k = nzb_w_inner(j,i)+2
1625	flux_d = flux_t(k-1)
1626	diss_d = diss_t(k-1)
1627	w_comp = w(k+1,j,i) + w(k,j,i)
1628	flux_t(k) = w_comp * ( &
1629	7.0 * ( w(k+1,j,i) + w(k,j,i) ) &
1630	- ( w(k+2,j,i) + w(k-1,j,i) ) ) * adv_mom_3
1631	diss_t(k) = - ABS( w_comp ) * ( &
1632	3.0 * ( w(k+1,j,i) - w(k,j,i) ) &
1633	- ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
1634
1635	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1636	- flux_d - diss_d ) * ddzu(k+1)
1637	!
1638	!-- WS5
1639	DO k = nzb_w_inner(j,i)+3, nzt-2
1640	flux_d = flux_t(k-1)
1641	diss_d = diss_t(k-1)
1642	w_comp = w(k+1,j,i) + w(k,j,i)
1643	flux_t(k) = w_comp * ( &
1644	37.0 * ( w(k+1,j,i) + w(k,j,i) ) &
1645	- 8.0 * ( w(k+2,j,i) + w(k-1,j,i) ) &
1646	+ ( w(k+3,j,i) + w(k-2,j,i) ) ) * adv_mom_5
1647	diss_t(k) = - ABS( w_comp ) * ( &
1648	10.0 * ( w(k+1,j,i) - w(k,j,i) ) &
1649	- 5.0 * ( w(k+2,j,i) - w(k-1,j,i) ) &
1650	+ ( w(k+3,j,i) - w(k-2,j,i) ) ) * adv_mom_5
1651
1652	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1653	- flux_d - diss_d ) * ddzu(k+1)
1654	ENDDO
1655	!-- WS3 as an intermediate step (top)
1656	k = nzt - 1
1657	flux_d = flux_t(k-1)
1658	diss_d = diss_t(k-1)
1659	w_comp = w(k+1,j,i) + w(k,j,i)
1660	flux_t(k) = w_comp * ( &
1661	7.0 * ( w(k+1,j,i) + w(k,j,i) ) &
1662	- ( w(k+2,j,i) + w(k-1,j,i) ) ) *adv_mom_3
1663	diss_t(k) = - ABS( w_comp ) * ( &
1664	3.0 * ( w(k+1,j,i) - w(k,j,i) ) &
1665	- ( w(k+2,j,i) - w(k-1,j,i) ) ) * adv_mom_3
1666
1667	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1668	- flux_d - diss_d ) * ddzu(k+1)
1669	!
1670	!-- 2nd order scheme (top)
1671	k = nzt
1672	flux_d = flux_t(k-1)
1673	diss_d = diss_t(k-1)
1674	w_comp = w(k+1,j,i) + w(k,j,i)
1675	flux_t(k) = w_comp * ( w(k+1,j,i) + w(k,j,i) ) * 0.25
1676	diss_t(k) = diss_2nd( w(k+1,j,i), w(k+1,j,i), w(k,j,i), w(k-1,j,i), &
1677	w(k-2,j,i), w_comp, 0.25, ddzu(k+1) )
1678
1679	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1680	- flux_d - diss_d ) * ddzu(k+1)
1681
1682	DO k = nzb_w_inner(j,i), nzt
1683	sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &
1684	+ ( flux_t(k) + diss_t(k) ) &
1685	* weight_substep(intermediate_timestep_count)
1686	ENDDO
1687
1688	END SUBROUTINE advec_w_ws_ij
1689
1690
1691	!------------------------------------------------------------------------------!
1692	! Scalar advection - Call for all grid points
1693	!------------------------------------------------------------------------------!
1694	SUBROUTINE advec_s_ws( sk, sk_char )
1695
1696	USE arrays_3d
1697	USE constants
1698	USE control_parameters
1699	USE grid_variables
1700	USE indices
1701	USE statistics
1702
1703	IMPLICIT NONE
1704
1705	INTEGER :: i, j, k, tn = 0
1706	REAL, DIMENSION(:,:,:), POINTER :: sk
1707	REAL :: flux_d, diss_d, u_comp, v_comp
1708	REAL, DIMENSION(nzb:nzt+1) :: flux_r, diss_r, flux_n, diss_n
1709	REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local, swap_diss_y_local, &
1710	flux_t, diss_t
1711	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local, &
1712	swap_diss_x_local
1713	CHARACTER (LEN = *), INTENT(IN) :: sk_char
1714
1715	!
1716	!-- Compute the fluxes for the whole left boundary of the processor domain.
1717	i = nxl
1718	DO j = nys, nyn
1719	IF ( boundary_flags(j,i) == 6 ) THEN
1720
1721	DO k = nzb_s_inner(j,i)+1, nzt
1722	u_comp = u(k,j,i) - u_gtrans
1723	swap_flux_x_local(k,j) = u_comp * ( &
1724	sk(k,j,i) + sk(k,j,i-1)) * 0.5
1725	swap_diss_x_local(k,j) = diss_2nd( sk(k,j,i+2), sk(k,j,i+1), &
1726	sk(k,j,i), sk(k,j,i-1), &
1727	sk(k,j,i-1), u_comp, &
1728	0.5, ddx )
1729	ENDDO
1730
1731	ELSE
1732
1733	DO k = nzb_s_inner(j,i)+1, nzt
1734	u_comp = u(k,j,i) - u_gtrans
1735	swap_flux_x_local(k,j) = u_comp*( &
1736	37.0 * ( sk(k,j,i)+sk(k,j,i-1) ) &
1737	- 8.0 * ( sk(k,j,i+1) + sk(k,j,i-2) ) &
1738	+ ( sk(k,j,i+2) + sk(k,j,i-3) ) )&
1739	* adv_sca_5
1740	swap_diss_x_local(k,j) = - ABS( u_comp ) * ( &
1741	10.0 * (sk(k,j,i) - sk(k,j,i-1) ) &
1742	- 5.0 * ( sk(k,j,i+1) - sk(k,j,i-2) ) &
1743	+ ( sk(k,j,i+2) - sk(k,j,i-3) ) )&
1744	* adv_sca_5
1745	ENDDO
1746	ENDIF
1747	ENDDO
1748	!
1749	!-- The following loop computes the horizontal fluxes for the interior of the
1750	!-- processor domain plus south boundary points. Furthermore tendency terms
1751	!-- are computed.
1752	DO i = nxl, nxr
1753	j = nys
1754	IF ( boundary_flags(j,i) == 8 ) THEN
1755
1756	DO k = nzb_s_inner(j,i)+1, nzt
1757	v_comp = v(k,j,i) - v_gtrans
1758	swap_flux_y_local(k) = v_comp * &
1759	( sk(k,j,i) + sk(k,j-1,i) ) * 0.5
1760	swap_diss_y_local(k) = diss_2nd( sk(k,j+2,i), sk(k,j+1,i), &
1761	sk(k,j,i), sk(k,j-1,i), &
1762	sk(k,j-1,i), v_comp, 0.5, ddy )
1763	ENDDO
1764
1765	ELSE
1766
1767	DO k = nzb_s_inner(j,i)+1, nzt
1768	v_comp = v(k,j,i) - v_gtrans
1769	swap_flux_y_local(k) = v_comp * ( &
1770	37.0 * ( sk(k,j,i) + sk(k,j-1,i) ) &
1771	- 8.0 * ( sk(k,j+1,i) + sk(k,j-2,i) ) &
1772	+ ( sk(k,j+2,i) + sk(k,j-3,i) ) ) &
1773	* adv_sca_5
1774	swap_diss_y_local(k)= - ABS( v_comp ) * ( &
1775	10.0 * ( sk(k,j,i) - sk(k,j-1,i) ) &
1776	- 5.0 * ( sk(k,j+1,i) - sk(k,j-2,i) ) &
1777	+ ( sk(k,j+2,i)-sk(k,j-3,i) ) ) &
1778	* adv_sca_5
1779	ENDDO
1780
1781	ENDIF
1782
1783	DO j = nys, nyn
1784	IF ( boundary_flags(j,i) /= 0 ) THEN
1785	!
1786	!-- Degrade the order for Dirichlet bc. at the outflow boundary
1787	SELECT CASE ( boundary_flags(j,i) )
1788
1789	CASE ( 1 )
1790	DO k = nzb_s_inner(j,i)+1, nzt
1791	u_comp = u(k,j,i+1) - u_gtrans
1792	flux_r(k) = u_comp * ( &
1793	7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
1794	- ( sk(k,j,i+2) + sk(k,j,i-1) ) ) &
1795	* adv_sca_3
1796	diss_r(k) = - ABS( u_comp ) * ( &
1797	3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
1798	- ( sk(k,j,i+2) - sk(k,j,i-1) ) ) &
1799	* adv_sca_3
1800	ENDDO
1801
1802	CASE ( 2 )
1803	DO k = nzb_s_inner(j,i)+1, nzt
1804	u_comp = u(k,j,i+1) - u_gtrans
1805	flux_r(k) = u_comp * ( sk(k,j,i+1) + sk(k,j,i) ) * 0.5
1806	diss_r(k) = diss_2nd( sk(k,j,i+1), sk(k,j,i+1), &
1807	sk(k,j,i), sk(k,j,i-1), &
1808	sk(k,j,i-2), u_comp, 0.5, ddx )
1809	ENDDO
1810
1811	CASE ( 3 )
1812	DO k = nzb_s_inner(j,i)+1, nzt
1813	v_comp = v(k,j+1,i) - v_gtrans
1814	flux_n(k) = v_comp * ( &
1815	7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
1816	- ( sk(k,j+2,i) + sk(k,j-1,i) ) ) &
1817	* adv_sca_3
1818	diss_n(k) = - ABS( v_comp ) * ( &
1819	3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
1820	- ( sk(k,j+2,i) - sk(k,j-1,i) ) ) &
1821	* adv_sca_3
1822	ENDDO
1823
1824	CASE ( 4 )
1825	DO k = nzb_s_inner(j,i)+1, nzt
1826	v_comp = v(k,j+1,i) - v_gtrans
1827	flux_n(k) = v_comp * ( sk(k,j+1,i) + sk(k,j,i) ) * 0.5
1828	diss_n(k) = diss_2nd( sk(k,j+1,i), sk(k,j+1,i), &
1829	sk(k,j,i), sk(k,j-1,i), &
1830	sk(k,j-2,i), v_comp, 0.5, ddy )
1831	ENDDO
1832
1833	CASE ( 5 )
1834	DO k = nzb_w_inner(j,i)+1, nzt
1835	u_comp = u(k,j,i+1) - u_gtrans
1836	flux_r(k) = u_comp * ( &
1837	7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
1838	- ( sk(k,j,i+2) + sk(k,j,i-1) ) ) &
1839	* adv_sca_3
1840	diss_r(k) = - ABS( u_comp ) * ( &
1841	3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
1842	- ( sk(k,j,i+2) - sk(k,j,i-1) ) ) &
1843	* adv_sca_3
1844	ENDDO
1845
1846	CASE ( 6 )
1847	DO k = nzb_s_inner(j,i)+1, nzt
1848	u_comp = u(k,j,i+1) - u_gtrans
1849	flux_r(k) = u_comp * ( &
1850	7.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
1851	- ( sk(k,j,i+2) + sk(k,j,i-1) ) ) &
1852	* adv_sca_3
1853	diss_r(k) = - ABS( u_comp ) * ( &
1854	3.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
1855	- ( sk(k,j,i+2) - sk(k,j,i-1) ) ) &
1856	* adv_sca_3
1857	ENDDO
1858
1859	CASE ( 7 )
1860	DO k = nzb_s_inner(j,i)+1, nzt
1861	v_comp = v(k,j+1,i) - v_gtrans
1862	flux_n(k) = v_comp * ( &
1863	7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
1864	- ( sk(k,j+2,i) + sk(k,j-1,i) ) ) &
1865	* adv_sca_3
1866	diss_n(k) = - ABS( v_comp ) * ( &
1867	3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
1868	- ( sk(k,j+2,i) - sk(k,j-1,i) ) ) &
1869	* adv_sca_3
1870	ENDDO
1871
1872	CASE ( 8 )
1873	DO k = nzb_s_inner(j,i)+1, nzt
1874	v_comp = v(k,j+1,i) - v_gtrans
1875	flux_n(k) = v_comp * ( &
1876	7.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
1877	- ( sk(k,j+2,i) + sk(k,j-1,i) ) ) &
1878	* adv_sca_3
1879	diss_n(k) = - ABS( v_comp ) * ( &
1880	3.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
1881	- ( sk(k,j+2,i) - sk(k,j-1,i) ) ) &
1882	* adv_sca_3
1883	ENDDO
1884
1885	CASE DEFAULT
1886
1887	END SELECT
1888
1889	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
1890	boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) &
1891	THEN
1892
1893	DO k = nzb_s_inner(j,i)+1, nzt
1894	v_comp = v(k,j+1,i) - v_gtrans
1895	flux_n(k) = v_comp * ( &
1896	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
1897	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
1898	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) &
1899	* adv_sca_5
1900	diss_n(k) = - ABS( v_comp ) * ( &
1901	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
1902	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
1903	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) &
1904	* adv_sca_5
1905	ENDDO
1906
1907	ELSE
1908
1909	DO k = nzb_s_inner(j,i)+1, nzt
1910	u_comp = u(k,j,i+1) - u_gtrans
1911	flux_r(k) = u_comp * ( &
1912	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
1913	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
1914	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) &
1915	* adv_sca_5
1916	diss_r(k) = - ABS( u_comp ) * ( &
1917	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
1918	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
1919	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) &
1920	* adv_sca_5
1921	ENDDO
1922
1923	ENDIF
1924
1925	ELSE
1926
1927	DO k = nzb_s_inner(j,i)+1, nzt
1928	u_comp = u(k,j,i+1) - u_gtrans
1929	flux_r(k) = u_comp * ( &
1930	37.0 * ( sk(k,j,i+1) + sk(k,j,i) ) &
1931	- 8.0 * ( sk(k,j,i+2) + sk(k,j,i-1) ) &
1932	+ ( sk(k,j,i+3) + sk(k,j,i-2) ) ) &
1933	* adv_sca_5
1934	diss_r(k) = - ABS( u_comp ) * ( &
1935	10.0 * ( sk(k,j,i+1) - sk(k,j,i) ) &
1936	- 5.0 * ( sk(k,j,i+2) - sk(k,j,i-1) ) &
1937	+ ( sk(k,j,i+3) - sk(k,j,i-2) ) ) &
1938	* adv_sca_5
1939
1940	v_comp = v(k,j+1,i) - v_gtrans
1941	flux_n(k) = v_comp * ( &
1942	37.0 * ( sk(k,j+1,i) + sk(k,j,i) ) &
1943	- 8.0 * ( sk(k,j+2,i) + sk(k,j-1,i) ) &
1944	+ ( sk(k,j+3,i) + sk(k,j-2,i) ) ) &
1945	* adv_sca_5
1946	diss_n(k) = - ABS( v_comp ) * ( &
1947	10.0 * ( sk(k,j+1,i) - sk(k,j,i) ) &
1948	- 5.0 * ( sk(k,j+2,i) - sk(k,j-1,i) ) &
1949	+ ( sk(k,j+3,i) - sk(k,j-2,i) ) ) &
1950	* adv_sca_5
1951	ENDDO
1952
1953	ENDIF
1954
1955	DO k = nzb_s_inner(j,i)+1, nzt
1956	tend(k,j,i) = tend(k,j,i) - ( &
1957	( flux_r(k) + diss_r(k) &
1958	- swap_flux_x_local(k,j) - swap_diss_x_local(k,j) ) * ddx &
1959	+ ( flux_n(k) + diss_n(k) &
1960	- swap_flux_y_local(k) - swap_diss_y_local(k) ) * ddy)
1961
1962	swap_flux_x_local(k,j) = flux_r(k)
1963	swap_diss_x_local(k,j) = diss_r(k)
1964	swap_flux_y_local(k) = flux_n(k)
1965	swap_diss_y_local(k) = diss_n(k)
1966	ENDDO
1967	ENDDO
1968	ENDDO
1969
1970	!
1971	!-- Vertical advection, degradation of order near surface and top.
1972	!-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
1973	!-- statistical evaluation the top flux at the surface should be 0
1974	DO i = nxl, nxr
1975	DO j = nys, nyn
1976	!
1977	!-- 2nd order scheme (bottom)
1978	k=nzb_s_inner(j,i)+1
1979	!
1980	!-- The fluxes flux_d and diss_d at the surface are 0. Due to static
1981	!-- reasons the top flux at the surface should be 0.
1982	flux_t(nzb_s_inner(j,i)) = 0.0
1983	diss_t(nzb_s_inner(j,i)) = 0.0
1984	flux_d = flux_t(k-1)
1985	diss_d = diss_t(k-1)
1986	flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
1987	diss_t(k) = diss_2nd( sk(k+2,j,i), sk(k+1,j,i), sk(k,j,i), &
1988	sk(k,j,i), sk(k,j,i), w(k,j,i), &
1989	0.5, ddzw(k) )
1990	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
1991	- flux_d - diss_d ) * ddzw(k)
1992	!
1993	!-- WS3 as an intermediate step (bottom)
1994	k = nzb_s_inner(j,i)+2
1995	flux_d = flux_t(k-1)
1996	diss_d = diss_t(k-1)
1997	flux_t(k) = w(k,j,i) * ( &
1998	7.0 * ( sk(k+1,j,i) + sk(k,j,i) ) &
1999	- ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
2000	diss_t(k) = - ABS( w(k,j,i) ) * ( &
2001	3.0 * ( sk(k+1,j,i) - sk(k,j,i) ) &
2002	- ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
2003	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2004	- flux_d - diss_d ) * ddzw(k)
2005	!
2006	!-- WS5
2007	DO k = nzb_s_inner(j,i)+3, nzt-2
2008	flux_d = flux_t(k-1)
2009	diss_d = diss_t(k-1)
2010	flux_t(k) = w(k,j,i) * ( &
2011	37.0 * ( sk(k+1,j,i) + sk(k,j,i) ) &
2012	- 8.0 * ( sk(k+2,j,i) + sk(k-1,j,i) ) &
2013	+ ( sk(k+3,j,i) + sk(k-2,j,i) ) ) * adv_sca_5
2014	diss_t(k) = - ABS(w(k,j,i)) * ( &
2015	10.0 * ( sk(k+1,j,i) -sk(k,j,i) ) &
2016	- 5.0 * ( sk(k+2,j,i) - sk(k-1,j,i) ) &
2017	+ ( sk(k+3,j,i) - sk(k-2,j,i) ) ) * adv_sca_5
2018
2019	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2020	- flux_d - diss_d ) * ddzw(k)
2021	ENDDO
2022	!
2023	!-- WS3 as an intermediate step (top)
2024	k = nzt - 1
2025	flux_d = flux_t(k-1)
2026	diss_d = diss_t(k-1)
2027	flux_t(k) = w(k,j,i) * ( &
2028	7.0 * ( sk(k+1,j,i) + sk(k,j,i) ) &
2029	- ( sk(k+2,j,i) + sk(k-1,j,i) ) ) * adv_sca_3
2030	diss_t(k) = - ABS(w(k,j,i)) * ( &
2031	3.0 * ( sk(k+1,j,i) - sk(k,j,i) ) &
2032	- ( sk(k+2,j,i) - sk(k-1,j,i) ) ) * adv_sca_3
2033
2034	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2035	- flux_d - diss_d ) * ddzw(k)
2036	!
2037	!-- 2nd order scheme (top)
2038	k = nzt
2039	flux_d = flux_t(k-1)
2040	diss_d = diss_t(k-1)
2041	flux_t(k) = w(k,j,i) * ( sk(k+1,j,i) + sk(k,j,i) ) * 0.5
2042	diss_t(k) = diss_2nd( sk(k+1,j,i), sk(k+1,j,i), sk(k,j,i), &
2043	sk(k-1,j,i), sk(k-2,j,i), w(k,j,i), &
2044	0.5, ddzw(k) )
2045
2046	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2047	- flux_d - diss_d ) * ddzw(k)
2048	!
2049	!-- evaluation of statistics
2050	SELECT CASE ( sk_char )
2051
2052	CASE ( 'pt' )
2053	DO k = nzb_s_inner(j,i), nzt
2054	sums_wspts_ws_l(k,tn) = sums_wspts_ws_l(k,tn) &
2055	+ ( flux_t(k) + diss_t(k) ) &
2056	* weight_substep(intermediate_timestep_count)
2057	ENDDO
2058	CASE ( 'sa' )
2059	DO k = nzb_s_inner(j,i), nzt
2060	sums_wssas_ws_l(k,tn) = sums_wssas_ws_l(k,tn) &
2061	+ ( flux_t(k) + diss_t(k) ) &
2062	* weight_substep(intermediate_timestep_count)
2063	ENDDO
2064	CASE ( 'q' )
2065	DO k = nzb_s_inner(j,i), nzt
2066	sums_wsqs_ws_l(k,tn) = sums_wsqs_ws_l(k,tn) &
2067	+ ( flux_t(k) + diss_t(k) ) &
2068	* weight_substep(intermediate_timestep_count)
2069	ENDDO
2070
2071	END SELECT
2072	ENDDO
2073	ENDDO
2074
2075
2076	END SUBROUTINE advec_s_ws
2077
2078
2079	!------------------------------------------------------------------------------!
2080	! Advection of u - Call for all grid points
2081	!------------------------------------------------------------------------------!
2082	SUBROUTINE advec_u_ws
2083
2084	USE arrays_3d
2085	USE constants
2086	USE control_parameters
2087	USE grid_variables
2088	USE indices
2089	USE statistics
2090
2091	IMPLICIT NONE
2092
2093	INTEGER :: i, j, k, tn = 0
2094	REAL :: gu, gv, flux_d, diss_d, v_comp, w_comp
2095	REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_u, swap_diss_y_local_u
2096	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_u, &
2097	swap_diss_x_local_u
2098	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_r, diss_r, flux_n, &
2099	diss_n, u_comp
2100
2101	gu = 2.0 * u_gtrans
2102	gv = 2.0 * v_gtrans
2103
2104	!
2105	!-- Compute the fluxes for the whole left boundary of the processor domain.
2106	i = nxlu
2107	DO j = nys, nyn
2108	IF( boundary_flags(j,i) == 5 ) THEN
2109
2110	DO k = nzb_u_inner(j,i)+1, nzt
2111	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
2112	swap_flux_x_local_u(k,j) = u_comp(k) * &
2113	( u(k,j,i) + u(k,j,i-1) ) * 0.25
2114	swap_diss_x_local_u(k,j) = diss_2nd( u(k,j,i+2), u(k,j,i+1), &
2115	u(k,j,i), u(k,j,i-1), &
2116	u(k,j,i-1), u_comp(k), &
2117	0.25, ddx )
2118	ENDDO
2119
2120	ELSE
2121
2122	DO k = nzb_u_inner(j,i)+1, nzt
2123	u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu
2124	swap_flux_x_local_u(k,j) = u_comp(k) * ( &
2125	37.0 * ( u(k,j,i) + u(k,j,i-1) ) &
2126	- 8.0 * ( u(k,j,i+1) + u(k,j,i-2) ) &
2127	+ (u(k,j,i+2)+u(k,j,i-3) ) ) &
2128	* adv_mom_5
2129	swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( &
2130	10.0 * ( u(k,j,i) - u(k,j,i-1) ) &
2131	- 5.0 * ( u(k,j,i+1) - u(k,j,i-2) ) &
2132	+ ( u(k,j,i+2) - u(k,j,i-3) ) )&
2133	* adv_mom_5
2134	ENDDO
2135
2136	ENDIF
2137
2138	ENDDO
2139
2140	DO i = nxlu, nxr
2141	!
2142	!-- The following loop computes the fluxes for the south boundary points
2143	j = nys
2144	IF ( boundary_flags(j,i) == 8 ) THEN
2145	!
2146	!-- Compute southside fluxes for the south boundary of PE domain
2147	DO k = nzb_u_inner(j,i)+1, nzt
2148	v_comp = v(k,j,i) + v(k,j,i-1) - gv
2149	swap_flux_y_local_u(k) = v_comp * &
2150	( u(k,j,i) + u(k,j-1,i) ) * 0.25
2151	swap_diss_y_local_u(k) = diss_2nd( u(k,j+2,i), u(k,j+1,i), &
2152	u(k,j,i), u(k,j-1,i), &
2153	u(k,j-1,i), v_comp, &
2154	0.25, ddy )
2155	ENDDO
2156
2157	ELSE
2158
2159	DO k = nzb_u_inner(j,i)+1, nzt
2160	v_comp = v(k,j,i) + v(k,j,i-1) - gv
2161	swap_flux_y_local_u(k) = v_comp * ( &
2162	37.0 * ( u(k,j,i) + u(k,j-1,i) ) &
2163	- 8.0 * ( u(k,j+1,i) + u(k,j-2,i) ) &
2164	+ ( u(k,j+2,i) + u(k,j-3,i) ) ) &
2165	* adv_mom_5
2166	swap_diss_y_local_u(k) = - ABS( v_comp ) * ( &
2167	10.0 * ( u(k,j,i) - u(k,j-1,i) ) &
2168	- 5.0 * ( u(k,j+1,i) - u(k,j-2,i) ) &
2169	+ ( u(k,j+2,i) - u(k,j-3,i) ) ) &
2170	* adv_mom_5
2171	ENDDO
2172
2173	ENDIF
2174	!
2175	!-- Computation of interior fluxes and tendency terms
2176	DO j = nys, nyn
2177	IF ( boundary_flags(j,i) /= 0 ) THEN
2178	!
2179	!-- Degrade the order for Dirichlet bc. at the outflow boundary
2180	SELECT CASE ( boundary_flags(j,i) )
2181
2182	CASE ( 1 )
2183	DO k = nzb_u_inner(j,i)+1, nzt
2184	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2185	flux_r(k) = ( u_comp(k) - gu ) * ( &
2186	7.0 * ( u(k,j,i+1) + u(k,j,i) ) &
2187	- ( u(k,j,i+2) + u(k,j,i-1) ) ) &
2188	* adv_mom_3
2189	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2190	3.0 * ( u(k,j,i+1) - u(k,j,i) ) &
2191	- ( u(k,j,i+2) - u(k,j,i-1) ) ) &
2192	* adv_mom_3
2193	ENDDO
2194
2195	CASE ( 2 )
2196	DO k = nzb_u_inner(j,i)+1, nzt
2197	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2198	flux_r(k) = ( u_comp(k) - gu ) * &
2199	( u(k,j,i+1) + u(k,j,i) ) * 0.25
2200	diss_r(k) = diss_2nd( u(k,j,i+1), u(k,j,i+1), &
2201	u(k,j,i), u(k,j,i-1), &
2202	u(k,j,i-2), u_comp(k) ,0.25 ,ddx)
2203	ENDDO
2204
2205	CASE ( 3 )
2206	DO k = nzb_u_inner(j,i)+1, nzt
2207	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2208	flux_n(k) = v_comp * ( &
2209	7.0 * ( u(k,j+1,i) + u(k,j,i) ) &
2210	- ( u(k,j+2,i) + u(k,j-1,i) ) ) &
2211	* adv_mom_3
2212	diss_n(k) = - ABS( v_comp ) * ( &
2213	3.0 * ( u(k,j+1,i) - u(k,j,i) ) &
2214	- ( u(k,j+2,i) - u(k,j-1,i) ) ) &
2215	* adv_mom_3
2216	ENDDO
2217
2218	CASE ( 4 )
2219	DO k = nzb_u_inner(j,i)+1, nzt
2220	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2221	flux_n(k) = v_comp * ( u(k,j+1,i) + u(k,j,i) ) * 0.25
2222	diss_n(k) = diss_2nd( u(k,j+1,i), u(k,j+1,i), &
2223	u(k,j,i), u(k,j-1,i), &
2224	u(k,j-2,i), v_comp, 0.25, ddy )
2225	ENDDO
2226
2227	CASE ( 5 )
2228	DO k = nzb_u_inner(j,i)+1, nzt
2229	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2230	flux_r(k) = ( u_comp(k) - gu ) * ( &
2231	7.0 * ( u(k,j,i+1) + u(k,j,i) ) &
2232	- ( u(k,j,i+2) + u(k,j,i-1) ) ) &
2233	* adv_mom_3
2234	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2235	3.0 * ( u(k,j,i+1) - u(k,j,i) ) &
2236	- ( u(k,j,i+2) - u(k,j,i-1) ) ) &
2237	* adv_mom_3
2238	ENDDO
2239
2240	CASE ( 7 )
2241	DO k = nzb_u_inner(j,i)+1, nzt
2242	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2243	flux_n(k) = v_comp * ( &
2244	7.0 * ( u(k,j+1,i) + u(k,j,i) ) &
2245	- ( u(k,j+2,i) + u(k,j-1,i) ) ) &
2246	* adv_mom_3
2247	diss_n(k) = - ABS( v_comp ) * ( &
2248	3.0 * ( u(k,j+1,i) - u(k,j,i) ) &
2249	- ( u(k,j+2,i) - u(k,j-1,i) ) ) &
2250	* adv_mom_3
2251	ENDDO
2252
2253	CASE ( 8 )
2254	DO k = nzb_u_inner(j,i)+1, nzt
2255	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2256	flux_n(k) = v_comp * ( &
2257	7.0 * ( u(k,j+1,i) + u(k,j,i) ) &
2258	- ( u(k,j+2,i) + u(k,j-1,i) ) ) &
2259	* adv_mom_3
2260	diss_n(k) = - ABS( v_comp ) * ( &
2261	3.0 * ( u(k,j+1,i) - u(k,j,i) ) &
2262	- ( u(k,j+2,i) - u(k,j-1,i) ) ) &
2263	* adv_mom_3
2264	ENDDO
2265
2266	CASE DEFAULT
2267
2268	END SELECT
2269	!
2270	!-- Compute the crosswise 5th order fluxes at the outflow
2271	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
2272	boundary_flags(j,i) == 5 ) THEN
2273
2274	DO k = nzb_u_inner(j,i)+1, nzt
2275	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2276	flux_n(k) = v_comp * ( &
2277	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
2278	- 8.0 * ( u(k,j+2,i) +u(k,j-1,i) ) &
2279	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) &
2280	* adv_mom_5
2281	diss_n(k) = - ABS( v_comp ) * ( &
2282	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
2283	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
2284	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) &
2285	* adv_mom_5
2286	ENDDO
2287
2288	ELSE
2289
2290	DO k = nzb_u_inner(j,i)+1, nzt
2291	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2292	flux_r(k) = ( u_comp(k) - gu ) * ( &
2293	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
2294	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
2295	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) &
2296	* adv_mom_5
2297	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2298	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
2299	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
2300	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) &
2301	* adv_mom_5
2302	ENDDO
2303
2304	ENDIF
2305
2306	ELSE
2307
2308	DO k = nzb_u_inner(j,i)+1, nzt
2309	u_comp(k) = u(k,j,i+1) + u(k,j,i)
2310	flux_r(k) = ( u_comp(k) - gu ) * ( &
2311	37.0 * ( u(k,j,i+1) + u(k,j,i) ) &
2312	- 8.0 * ( u(k,j,i+2) + u(k,j,i-1) ) &
2313	+ ( u(k,j,i+3) + u(k,j,i-2) ) ) &
2314	* adv_mom_5
2315	diss_r(k) = - ABS( u_comp(k) - gu ) * ( &
2316	10.0 * ( u(k,j,i+1) - u(k,j,i) ) &
2317	- 5.0 * ( u(k,j,i+2) - u(k,j,i-1) ) &
2318	+ ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5
2319
2320	v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv
2321	flux_n(k) = v_comp * ( &
2322	37.0 * ( u(k,j+1,i) + u(k,j,i) ) &
2323	- 8.0 * ( u(k,j+2,i) + u(k,j-1,i) ) &
2324	+ ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5
2325	diss_n(k) = - ABS( v_comp ) * ( &
2326	10.0 * ( u(k,j+1,i) - u(k,j,i) ) &
2327	- 5.0 * ( u(k,j+2,i) - u(k,j-1,i) ) &
2328	+ ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5
2329
2330	ENDDO
2331
2332	ENDIF
2333
2334	DO k = nzb_u_inner(j,i)+1, nzt
2335
2336	tend(k,j,i) = tend(k,j,i) - ( &
2337	( flux_r(k) + diss_r(k) &
2338	- swap_flux_x_local_u(k,j) - swap_diss_x_local_u(k,j) ) * ddx &
2339	+ ( flux_n(k) + diss_n(k) &
2340	- swap_flux_y_local_u(k) - swap_diss_y_local_u(k) ) * ddy )
2341
2342	swap_flux_x_local_u(k,j) = flux_r(k)
2343	swap_diss_x_local_u(k,j) = diss_r(k)
2344	swap_flux_y_local_u(k) = flux_n(k)
2345	swap_diss_y_local_u(k) = diss_n(k)
2346
2347	sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &
2348	+ ( flux_r(k) &
2349	* ( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2350	/ ( u_comp(k) - gu + 1.0E-20 ) &
2351	+ diss_r(k) &
2352	* ABS( u_comp(k) - 2.0 * hom(k,1,1,0) ) &
2353	/ ( ABS( u_comp(k) - gu) + 1.0E-20) ) &
2354	* weight_substep(intermediate_timestep_count)
2355	ENDDO
2356	sums_us2_ws_l(nzb_u_inner(j,i),tn) = sums_us2_ws_l(nzb_u_inner(j,i)+1,tn)
2357	ENDDO
2358	ENDDO
2359
2360	!
2361	!-- Vertical advection, degradation of order near surface and top.
2362	!-- The fluxes flux_d and diss_d at the surface are 0. Due to reasons of
2363	!-- statistical evaluation the top flux at the surface should be 0
2364	DO i = nxlu, nxr
2365	DO j = nys, nyn
2366	k = nzb_u_inner(j,i)+1
2367	!
2368	!-- The fluxes flux_d and diss_d at the surface are 0. Due to static
2369	!-- reasons the top flux at the surface should be 0.
2370	flux_t(nzb_u_inner(j,i)) = 0.0
2371	diss_t(nzb_u_inner(j,i)) = 0.0
2372	flux_d = flux_t(k-1)
2373	diss_d = diss_t(k-1)
2374	!
2375	!-- 2nd order scheme (bottom)
2376	w_comp = w(k,j,i) + w(k,j,i-1)
2377	flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
2378	diss_t(k) = diss_2nd( u(k+2,j,i), u(k+1,j,i), u(k,j,i), &
2379	0.0, 0.0, w_comp, 0.25, ddzw(k) )
2380
2381	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2382	- flux_d - diss_d ) * ddzw(k)
2383	!
2384	!-- WS3 as an intermediate step (bottom)
2385	k = nzb_u_inner(j,i)+2
2386	flux_d = flux_t(k-1)
2387	diss_d = diss_t(k-1)
2388	w_comp = w(k,j,i) + w(k,j,i-1)
2389	flux_t(k) = w_comp * ( &
2390	7.0 * ( u(k+1,j,i) + u(k,j,i) ) &
2391	- ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
2392	diss_t(k) = - ABS( w_comp ) * ( &
2393	3.0 * ( u(k+1,j,i) - u(k,j,i) ) &
2394	- ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
2395
2396	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2397	- flux_d - diss_d ) * ddzw(k)
2398	!
2399	!WS5
2400	DO k = nzb_u_inner(j,i)+3, nzt-2
2401
2402	flux_d = flux_t(k-1)
2403	diss_d = diss_t(k-1)
2404	w_comp = w(k,j,i) + w(k,j,i-1)
2405	flux_t(k) = w_comp * ( &
2406	37.0 * ( u(k+1,j,i) + u(k,j,i) ) &
2407	- 8.0 * ( u(k+2,j,i) + u(k-1,j,i) ) &
2408	+ ( u(k+3,j,i) + u(k-2,j,i) ) ) * adv_mom_5
2409	diss_t(k) = - ABS( w_comp ) * ( &
2410	10.0 * ( u(k+1,j,i) - u(k,j,i) ) &
2411	- 5.0 * ( u(k+2,j,i) - u(k-1,j,i) ) &
2412	+ ( u(k+3,j,i) - u(k-2,j,i) ) ) * adv_mom_5
2413
2414	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2415	- flux_d - diss_d ) * ddzw(k)
2416
2417	ENDDO
2418	!
2419	!-- WS3 as an intermediate step (top)
2420	k = nzt-1
2421	flux_d = flux_t(k-1)
2422	diss_d = diss_t(k-1)
2423	w_comp = w(k,j,i) + w(k,j,i-1)
2424	flux_t(k) = w_comp * ( &
2425	7.0 * ( u(k+1,j,i) + u(k,j,i) ) &
2426	- ( u(k+2,j,i) + u(k-1,j,i) ) ) * adv_mom_3
2427	diss_t(k) = - ABS( w_comp ) * ( &
2428	3.0 * ( u(k+1,j,i) - u(k,j,i) ) &
2429	- ( u(k+2,j,i) - u(k-1,j,i) ) ) * adv_mom_3
2430
2431	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2432	- flux_d - diss_d ) * ddzw(k)
2433	!
2434	!-- 2nd order scheme (top)
2435	k = nzt
2436	flux_d = flux_t(k-1)
2437	diss_d = diss_t(k-1)
2438	w_comp = w(k,j,i) + w(k,j,i-1)
2439	flux_t(k) = w_comp * ( u(k+1,j,i) + u(k,j,i) ) * 0.25
2440	diss_t(k) = diss_2nd( u(nzt+1,j,i), u(nzt+1,j,i), u(k,j,i), &
2441	u(k-1,j,i), u(k-2,j,i), w_comp, &
2442	0.25, ddzw(k))
2443
2444	tend(k,j,i) = tend(k,j,i) - ( flux_t(k) + diss_t(k) &
2445	- flux_d - diss_d ) * ddzw(k)
2446	!
2447	!-- at last vertical momentum flux is accumulated
2448	DO k = nzb_u_inner(j,i), nzt
2449	sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &
2450	+ ( flux_t(k) + diss_t(k) ) &
2451	* weight_substep(intermediate_timestep_count)
2452	ENDDO
2453	ENDDO
2454	ENDDO
2455
2456
2457	END SUBROUTINE advec_u_ws
2458
2459
2460	!------------------------------------------------------------------------------!
2461	! Advection of v - Call for all grid points
2462	!------------------------------------------------------------------------------!
2463	SUBROUTINE advec_v_ws
2464
2465	USE arrays_3d
2466	USE constants
2467	USE control_parameters
2468	USE grid_variables
2469	USE indices
2470	USE statistics
2471
2472	IMPLICIT NONE
2473
2474
2475	INTEGER :: i, j, k, tn = 0
2476	REAL :: gu, gv, flux_l, flux_s, flux_d, diss_l, diss_s, diss_d, &
2477	u_comp, w_comp
2478	REAL, DIMENSION(nzb+1:nzt) :: swap_flux_y_local_v, swap_diss_y_local_v
2479	REAL, DIMENSION(nzb+1:nzt,nys:nyn) :: swap_flux_x_local_v, &
2480	swap_diss_x_local_v
2481	REAL, DIMENSION(nzb:nzt+1) :: flux_t, diss_t, flux_n, diss_n, flux_r, &
2482	diss_r, v_comp
2483
2484	gu = 2.0 * u_gtrans
2485	gv = 2.0 * v_gtrans
2486	!
2487	!-- First compute the whole left boundary of the processor domain
2488	i = nxl
2489	DO j = nysv, nyn
2490
2491	IF ( boundary_flags(j,i) == 6 ) THEN
2492	DO k = nzb_v_inner(j,i)+1, nzt
2493	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2494	swap_flux_x_local_v(k,j) = u_comp * &
2495	( v(k,j,i) + v(k,j,i-1)) * 0.25
2496	swap_diss_x_local_v(k,j) = diss_2nd( v(k,j,i+2), v(k,j,i+1), &
2497	v(k,j,i), v(k,j,i-1), &
2498	v(k,j,i-1), u_comp, &
2499	0.25, ddx )
2500	ENDDO
2501
2502	ELSE
2503
2504	DO k = nzb_v_inner(j,i)+1, nzt
2505	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2506	swap_flux_x_local_v(k,j) = u_comp * ( &
2507	37.0 * ( v(k,j,i) + v(k,j,i-1) ) &
2508	- 8.0 * ( v(k,j,i+1) + v(k,j,i-2) ) &
2509	+ ( v(k,j,i+2) + v(k,j,i-3) ) )&
2510	* adv_mom_5
2511	swap_diss_x_local_v(k,j) = - ABS( u_comp ) * ( &
2512	10.0 * ( v(k,j,i) - v(k,j,i-1) ) &
2513	- 5.0 * ( v(k,j,i+1) - v(k,j,i-2) ) &
2514	+ ( v(k,j,i+2) - v(k,j,i-3) ) )&
2515	* adv_mom_5
2516	ENDDO
2517
2518	ENDIF
2519
2520	ENDDO
2521
2522	DO i = nxl, nxr
2523	j = nysv
2524	IF ( boundary_flags(j,i) == 7 ) THEN
2525
2526	DO k = nzb_v_inner(j,i)+1, nzt
2527	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
2528	swap_flux_y_local_v(k) = v_comp(k) * &
2529	( v(k,j,i) + v(k,j-1,i) ) * 0.25
2530	swap_diss_y_local_v(k) = diss_2nd( v(k,j+2,i), v(k,j+1,i), &
2531	v(k,j,i), v(k,j-1,i), &
2532	v(k,j-1,i), v_comp(k), &
2533	0.25, ddy )
2534	ENDDO
2535
2536	ELSE
2537
2538	DO k = nzb_v_inner(j,i)+1, nzt
2539	v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv
2540	swap_flux_y_local_v(k) = v_comp(k) * ( &
2541	37.0 * ( v(k,j,i) + v(k,j-1,i) ) &
2542	- 8.0 * ( v(k,j+1,i) + v(k,j-2,i) ) &
2543	+ ( v(k,j+2,i) + v(k,j-3,i) ) ) &
2544	* adv_mom_5
2545	swap_diss_y_local_v(k) = - ABS( v_comp(k) ) * ( &
2546	10.0 * ( v(k,j,i) - v(k,j-1,i) ) &
2547	- 5.0 * ( v(k,j+1,i) - v(k,j-2,i) ) &
2548	+ ( v(k,j+2,i) - v(k,j-3,i) ) ) &
2549	* adv_mom_5
2550	ENDDO
2551
2552	ENDIF
2553
2554	DO j = nysv, nyn
2555	IF ( boundary_flags(j,i) /= 0 ) THEN
2556	!
2557	!-- Degrade the order for Dirichlet bc. at the outflow boundary
2558	SELECT CASE ( boundary_flags(j,i) )
2559
2560	CASE ( 1 )
2561	DO k = nzb_v_inner(j,i)+1, nzt
2562	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2563	flux_r(k) = u_comp * ( &
2564	7.0 * (v(k,j,i+1) + v(k,j,i) ) &
2565	- ( v(k,j,i+2) + v(k,j,i-1) ) ) &
2566	* adv_mom_3
2567	diss_r(k) = - ABS( u_comp ) * ( &
2568	3.0 * ( v(k,j,i+1) - v(k,j,i) ) &
2569	- ( v(k,j,i+2) - v(k,j,i-1) ) ) &
2570	* adv_mom_3
2571	ENDDO
2572
2573	CASE ( 2 )
2574	DO k = nzb_v_inner(j,i)+1, nzt
2575	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2576	flux_r(k) = u_comp * ( v(k,j,i+1) + v(k,j,i) ) * 0.25
2577	diss_r(k) = diss_2nd( v(k,j,i+1), v(k,j,i+1), &
2578	v(k,j,i), v(k,j,i-1), &
2579	v(k,j,i-2), u_comp, 0.25, ddx )
2580	ENDDO
2581
2582	CASE ( 3 )
2583	DO k = nzb_v_inner(j,i)+1, nzt
2584	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2585	flux_n(k) = ( v_comp(k)- gv ) * ( &
2586	7.0 * ( v(k,j+1,i) + v(k,j,i) ) &
2587	- ( v(k,j+2,i) + v(k,j-1,i) ) ) &
2588	* adv_mom_3
2589	diss_n(k) = - ABS(v_comp(k) - gv) * ( &
2590	3.0 * ( v(k,j+1,i) - v(k,j,i) ) &
2591	- ( v(k,j+2,i) - v(k,j-1,i) ) ) &
2592	* adv_mom_3
2593	ENDDO
2594
2595	CASE ( 4 )
2596	DO k = nzb_v_inner(j,i)+1, nzt
2597	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2598	flux_n(k) = ( v_comp(k) - gv ) * &
2599	( v(k,j+1,i) + v(k,j,i) ) * 0.25
2600	diss_n(k) = diss_2nd( v(k,j+1,i), v(k,j+1,i), &
2601	v(k,j,i), v(k,j-1,i), &
2602	v(k,j-2,i), v_comp(k), 0.25, ddy)
2603	ENDDO
2604
2605	CASE ( 5 )
2606	DO k = nzb_w_inner(j,i)+1, nzt
2607	u_comp = u(k,j-1,i) + u(k,j,i) - gu
2608	flux_r(k) = u_comp * ( &
2609	7.0 * ( v(k,j,i+1) + v(k,j,i) ) &
2610	- ( v(k,j,i+2) + v(k,j,i-1) ) ) &
2611	* adv_mom_3
2612	diss_r(k) = - ABS (u_comp ) * ( &
2613	3.0 * ( w(k,j,i+1) - w(k,j,i) ) &
2614	- ( v(k,j,i+2) - v(k,j,i-1) ) ) &
2615	* adv_mom_3
2616	ENDDO
2617
2618	CASE ( 6 )
2619	DO k = nzb_v_inner(j,i)+1, nzt
2620
2621	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2622	flux_r(k) = u_comp * ( &
2623	7.0 * ( v(k,j,i+1) + v(k,j,i) ) &
2624	- ( v(k,j,i+2) + v(k,j,i-1) ) ) &
2625	* adv_mom_3
2626	diss_r(k) = - ABS( u_comp ) * ( &
2627	3.0 * ( v(k,j,i+1) - v(k,j,i) ) &
2628	- ( v(k,j,i+2) - v(k,j,i-1) ) ) &
2629	* adv_mom_3
2630	ENDDO
2631
2632	CASE ( 7 )
2633	DO k = nzb_v_inner(j,i)+1, nzt
2634	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2635	flux_n(k) = ( v_comp(k) - gv ) * ( &
2636	7.0 * ( v(k,j+1,i) + v(k,j,i) ) &
2637	- ( v(k,j+2,i) + v(k,j-1,i) ) ) &
2638	* adv_mom_3
2639	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2640	3.0 * ( v(k,j+1,i) - v(k,j,i) ) &
2641	- ( v(k,j+2,i) - v(k,j-1,i) ) ) &
2642	* adv_mom_3
2643	ENDDO
2644
2645	CASE DEFAULT
2646
2647	END SELECT
2648	!
2649	!-- Compute the crosswise 5th order fluxes at the outflow
2650	IF ( boundary_flags(j,i) == 1 .OR. boundary_flags(j,i) == 2 .OR.&
2651	boundary_flags(j,i) == 5 .OR. boundary_flags(j,i) == 6 ) &
2652	THEN
2653
2654	DO k = nzb_v_inner(j,i)+1, nzt
2655	v_comp(k) = v(k,j+1,i) + v(k,j,i)
2656	flux_n(k) = ( v_comp(k) - gv ) * ( &
2657	37.0 * ( v(k,j+1,i) + v(k,j,i) ) &
2658	- 8.0 * ( v(k,j+2,i) + v(k,j-1,i) ) &
2659	+ ( v(k,j+3,i) + v(k,j-2,i) ) ) &
2660	* adv_mom_5
2661	diss_n(k) = - ABS( v_comp(k) - gv ) * ( &
2662	10.0 * ( v(k,j+1,i) - v(k,j,i) ) &
2663	- 5.0 * ( v(k,j+2,i) - v(k,j-1,i) ) &
2664	+ ( v(k,j+3,i) - v(k,j-2,i) ) ) &
2665	* adv_mom_5
2666	ENDDO
2667
2668	ELSE
2669
2670	DO k = nzb_v_inner(j,i)+1, nzt
2671	u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu
2672	flux_r(k) = u_comp * ( &
2673	37.0 * ( v(k,j,i+1) + v(k,j,i) ) &
2674	- 8.0 * ( v(k,j,i+2) + v(k,j,i-1) ) &
2675	+ ( v(k,j,i+3) + v(k,j,i-2)