Home

Context Navigation

advec_ws.f90 @ 711

Last change on this file since 711 was 711, checked in by raasch, 12 years ago
propset ID for advec_ws.f90
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 152.9 KB

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