Home

Context Navigation

source: palm/trunk/SOURCE/advec_ws.f90 @ 773

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

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |