Home

Context Navigation

advec_ws.f90 @ 714

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

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