Home

Context Navigation

advec_ws.f90 @ 802

Last change on this file since 802 was 802, checked in by suehring, 13 years ago
Last commit documented.
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 161.1 KB

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