Home

Context Navigation

advec_ws.f90 @ 706

Last change on this file since 706 was 706, checked in by suehring, 12 years ago
Last commit documented.
Property svn:executable set to ``*
File size: 151.8 KB

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