Home

Context Navigation

advec_ws.f90 @ 709

Last change on this file since 709 was 709, checked in by raasch, 12 years ago
formatting adjustments
Property svn:executable set to ``*
File size: 152.9 KB

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