Home

Context Navigation

advec_ws.f90 @ 673

Last change on this file since 673 was 673, checked in by suehring, 12 years ago
Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.
Property svn:executable set to ``*
File size: 151.4 KB

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