Home

Context Navigation

advec_ws.f90 @ 674

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

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