Home

Context Navigation

advec_ws.f90 @ 713

Last change on this file since 713 was 713, checked in by suehring, 12 years ago
reformatted advec_ws.f90, reformulate constants as broken numbers, bugfix in vertical advection of vertical velocity (concerning vector optimized routine)
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 159.9 KB

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