Home

Context Navigation

source: palm/trunk/SOURCE/poismg_fast_mod.f90 @ 1850

Last change on this file since 1850 was 1850, checked in by maronga, 8 years ago
added _mod string to several filenames to meet the naming convection for modules
Property svn:keywords set to `Id`
File size: 115.8 KB

Line
1	!> @file poismg_fast_mod.f90
2	!--------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the terms
6	! of the GNU General Public License as published by the Free Software Foundation,
7	! either version 3 of the License, or (at your option) any later version.
8	!
9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
12	!
13	! You should have received a copy of the GNU General Public License along with
14	! PALM. If not, see <http://www.gnu.org/licenses/>.
15	!
16	! Copyright 1997-2016 Leibniz Universitaet Hannover
17	!--------------------------------------------------------------------------------!
18	!
19	! Current revisions:
20	! -----------------
21	! Module renamed
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: poismg_fast_mod.f90 1850 2016-04-08 13:29:27Z maronga $
27	!
28	! 1762 2016-02-25 12:31:13Z hellstea
29	! Introduction of nested domain feature
30	!
31	! 1682 2015-10-07 23:56:08Z knoop
32	! Code annotations made doxygen readable
33	!
34	! 1609 2015-07-03 15:37:58Z maronga
35	! Bugfix: allow compilation without __parallel.
36	!
37	! 1575 2015-03-27 09:56:27Z raasch
38	! Initial revision.
39	! Routine re-written and optimised based on poismg.
40	!
41	! Following optimisations have been made:
42	! - vectorisation (for Intel-CPUs) of the red-black algorithm by resorting
43	! array elements with even and odd indices
44	! - explicit boundary conditions for building walls removed (solver is
45	! running through the buildings
46	! - reduced data transfer in case of ghost point exchange, because only
47	! "red" or "black" data points need to be exchanged. This is not applied
48	! for coarser grid levels, since for then the transfer time is latency bound
49	!
50	!
51	! Description:
52	! ------------
53	!> Solves the Poisson equation for the perturbation pressure with a multigrid
54	!> V- or W-Cycle scheme.
55	!>
56	!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
57	!> September 2000 - July 2001. It has been optimised for speed by Klaus
58	!> Ketelsen in November 2014.
59	!>
60	!> @attention Loop unrolling and cache optimization in SOR-Red/Black method
61	!> still does not give the expected speedup!
62	!>
63	!> @todo Further work required.
64	!------------------------------------------------------------------------------!
65	MODULE poismg_mod
66
67
68	USE cpulog, &
69	ONLY: cpu_log, log_point_s
70
71	USE kinds
72
73	USE pegrid
74
75	PRIVATE
76
77	INTEGER, SAVE :: ind_even_odd !< border index between even and odd k index
78	INTEGER, DIMENSION(:), SAVE, ALLOCATABLE :: even_odd_level !< stores ind_even_odd for all MG levels
79
80	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: f1_mg_b, f2_mg_b, f3_mg_b !< blocked version of f1_mg ...
81
82	INTERFACE poismg_fast
83	MODULE PROCEDURE poismg_fast
84	END INTERFACE poismg_fast
85
86	INTERFACE sort_k_to_even_odd_blocks
87	MODULE PROCEDURE sort_k_to_even_odd_blocks
88	MODULE PROCEDURE sort_k_to_even_odd_blocks_int
89	MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
90	END INTERFACE sort_k_to_even_odd_blocks
91
92	PUBLIC poismg_fast
93
94	CONTAINS
95
96	!------------------------------------------------------------------------------!
97	! Description:
98	! ------------
99	!> Solves the Poisson equation for the perturbation pressure with a multigrid
100	!> V- or W-Cycle scheme.
101	!------------------------------------------------------------------------------!
102	SUBROUTINE poismg_fast( r )
103
104	USE arrays_3d, &
105	ONLY: d, p_loc
106
107	USE control_parameters, &
108	ONLY: gathered_size, grid_level, grid_level_count, &
109	maximum_grid_level, message_string, mgcycles, mg_cycles, &
110	mg_switch_to_pe0_level, residual_limit, subdomain_size
111
112	USE cpulog, &
113	ONLY: cpu_log, log_point_s
114
115	USE indices, &
116	ONLY: nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,&
117	nyng, nyn_mg, nzb, nzt, nzt_mg
118
119	IMPLICIT NONE
120
121	REAL(wp) :: maxerror !<
122	REAL(wp) :: maximum_mgcycles !<
123	REAL(wp) :: residual_norm !<
124
125	REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !<
126
127	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !<
128
129
130	CALL cpu_log( log_point_s(29), 'poismg_fast', 'start' )
131	!
132	!-- Initialize arrays and variables used in this subroutine
133
134	!-- If the number of grid points of the gathered grid, which is collected
135	!-- on PE0, is larger than the number of grid points of an PE, than array
136	!-- p3 will be enlarged.
137	IF ( gathered_size > subdomain_size ) THEN
138	ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
139	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
140	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
141	mg_switch_to_pe0_level)+1) )
142	ELSE
143	ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
144	ENDIF
145
146	p3 = 0.0_wp
147
148
149	!
150	!-- Ghost boundaries have to be added to divergence array.
151	!-- Exchange routine needs to know the grid level!
152	grid_level = maximum_grid_level
153	CALL exchange_horiz( d, 1)
154	d(nzb,:,:) = d(nzb+1,:,:)
155
156	!
157	!-- Initiation of the multigrid scheme. Does n cycles until the
158	!-- residual is smaller than the given limit. The accuracy of the solution
159	!-- of the poisson equation will increase with the number of cycles.
160	!-- If the number of cycles is preset by the user, this number will be
161	!-- carried out regardless of the accuracy.
162	grid_level_count = 0
163	mgcycles = 0
164	IF ( mg_cycles == -1 ) THEN
165	maximum_mgcycles = 0
166	residual_norm = 1.0_wp
167	ELSE
168	maximum_mgcycles = mg_cycles
169	residual_norm = 0.0_wp
170	ENDIF
171
172	!
173	!-- Initial settings for sorting k-dimension from sequential order (alternate
174	!-- even/odd) into blocks of even and odd or vice versa
175	CALL init_even_odd_blocks
176
177	!
178	!-- Sort input arrays in even/odd blocks along k-dimension
179	CALL sort_k_to_even_odd_blocks( d, grid_level )
180	CALL sort_k_to_even_odd_blocks( p_loc, grid_level )
181
182	!
183	!-- The complete multigrid cycles are running in block mode, i.e. over
184	!-- seperate data blocks of even and odd indices
185	DO WHILE ( residual_norm > residual_limit .OR. &
186	mgcycles < maximum_mgcycles )
187
188	CALL next_mg_level_fast( d, p_loc, p3, r)
189
190	!
191	!-- Calculate the residual if the user has not preset the number of
192	!-- cycles to be performed
193	IF ( maximum_mgcycles == 0 ) THEN
194	CALL resid_fast( d, p_loc, r )
195	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
196
197	#if defined( __parallel )
198	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
199	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, &
200	MPI_SUM, comm2d, ierr)
201	#else
202	residual_norm = maxerror
203	#endif
204	residual_norm = SQRT( residual_norm )
205	ENDIF
206
207	mgcycles = mgcycles + 1
208
209	!
210	!-- If the user has not limited the number of cycles, stop the run in case
211	!-- of insufficient convergence
212	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
213	message_string = 'no sufficient convergence within 1000 cycles'
214	CALL message( 'poismg_fast', 'PA0283', 1, 2, 0, 6, 0 )
215	ENDIF
216
217	ENDDO
218
219	DEALLOCATE( p3 )
220	!
221	!-- Result has to be sorted back from even/odd blocks to sequential order
222	CALL sort_k_to_sequential( p_loc )
223	!
224	!-- Unset the grid level. Variable is used to determine the MPI datatypes for
225	!-- ghost point exchange
226	grid_level = 0
227
228	CALL cpu_log( log_point_s(29), 'poismg_fast', 'stop' )
229
230	END SUBROUTINE poismg_fast
231
232
233	!------------------------------------------------------------------------------!
234	! Description:
235	! ------------
236	!> Computes the residual of the perturbation pressure.
237	!------------------------------------------------------------------------------!
238	SUBROUTINE resid_fast( f_mg, p_mg, r )
239
240
241	USE arrays_3d, &
242	ONLY: f1_mg, f2_mg, f3_mg
243
244	USE control_parameters, &
245	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
246	inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l, &
247	nest_bound_n, nest_bound_r, nest_bound_s, outflow_l, &
248	outflow_n, outflow_r, outflow_s, topography
249
250	USE grid_variables, &
251	ONLY: ddx2_mg, ddy2_mg
252
253	USE indices, &
254	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, &
255	wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, &
256	wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, &
257	nys_mg, nyn_mg, nzb, nzt_mg
258
259	IMPLICIT NONE
260
261	INTEGER(iwp) :: i
262	INTEGER(iwp) :: j
263	INTEGER(iwp) :: k
264	INTEGER(iwp) :: l
265	INTEGER(iwp) :: km1 !<
266	INTEGER(iwp) :: kp1 !<
267
268	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
269	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
270	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
271	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
272	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
273	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
274	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
275	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
276	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
277
278	!
279	!-- Calculate the residual
280	l = grid_level
281
282	CALL cpu_log( log_point_s(53), 'resid', 'start' )
283	IF ( topography == 'flat' .OR. masking_method ) THEN
284	!$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
285	!$OMP DO
286	DO i = nxl_mg(l), nxr_mg(l)
287	DO j = nys_mg(l), nyn_mg(l)
288	!DIR$ IVDEP
289	DO k = ind_even_odd+1, nzt_mg(l)
290	km1 = k-ind_even_odd-1
291	kp1 = k-ind_even_odd
292	r(k,j,i) = f_mg(k,j,i) &
293	- ddx2_mg(l) * &
294	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
295	- ddy2_mg(l) * &
296	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
297	- f2_mg_b(k,l) * p_mg(kp1,j,i) &
298	- f3_mg_b(k,l) * p_mg(km1,j,i) &
299	+ f1_mg_b(k,l) * p_mg(k,j,i)
300	ENDDO
301	!DIR$ IVDEP
302	DO k = nzb+1, ind_even_odd
303	km1 = k+ind_even_odd
304	kp1 = k+ind_even_odd+1
305	r(k,j,i) = f_mg(k,j,i) &
306	- ddx2_mg(l) * &
307	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
308	- ddy2_mg(l) * &
309	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
310	- f2_mg_b(k,l) * p_mg(kp1,j,i) &
311	- f3_mg_b(k,l) * p_mg(km1,j,i) &
312	+ f1_mg_b(k,l) * p_mg(k,j,i)
313	ENDDO
314	ENDDO
315	ENDDO
316	!$OMP END PARALLEL
317	ELSE
318	!
319	!-- Choose flag array of this level
320	SELECT CASE ( l )
321	CASE ( 1 )
322	flags => wall_flags_1
323	CASE ( 2 )
324	flags => wall_flags_2
325	CASE ( 3 )
326	flags => wall_flags_3
327	CASE ( 4 )
328	flags => wall_flags_4
329	CASE ( 5 )
330	flags => wall_flags_5
331	CASE ( 6 )
332	flags => wall_flags_6
333	CASE ( 7 )
334	flags => wall_flags_7
335	CASE ( 8 )
336	flags => wall_flags_8
337	CASE ( 9 )
338	flags => wall_flags_9
339	CASE ( 10 )
340	flags => wall_flags_10
341	END SELECT
342
343	!$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
344	!$OMP DO
345	DO i = nxl_mg(l), nxr_mg(l)
346	DO j = nys_mg(l), nyn_mg(l)
347
348	!DIR$ IVDEP
349	DO k = ind_even_odd+1, nzt_mg(l)
350	km1 = k-ind_even_odd-1
351	kp1 = k-ind_even_odd
352
353	r(k,j,i) = f_mg(k,j,i) &
354	- ddx2_mg(l) * &
355	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 )) &
356	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
357	- ddy2_mg(l) &
358	* ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
359	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
360	- f2_mg(k,l) * p_mg(kp1,j,i) &
361	- f3_mg(k,l) * &
362	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 )) &
363	+ f1_mg(k,l) * p_mg(k,j,i)
364	ENDDO
365
366	!DIR$ IVDEP
367	DO k = nzb+1, ind_even_odd
368	km1 = k+ind_even_odd
369	kp1 = k+ind_even_odd+1
370
371	r(k,j,i) = f_mg(k,j,i) &
372	- ddx2_mg(l) * &
373	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 )) &
374	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
375	- ddy2_mg(l) &
376	* ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
377	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
378	- f2_mg(k,l) * p_mg(kp1,j,i) &
379	- f3_mg(k,l) * &
380	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 )) &
381	+ f1_mg(k,l) * p_mg(k,j,i)
382	ENDDO
383	!
384	!-- The residual within topography should be zero
385	WHERE( BTEST(flags(nzb+1:nzt_mg(l),j,i), 6) )
386	r(nzb+1:nzt_mg(l),j,i) = 0.0
387	END WHERE
388
389	ENDDO
390	ENDDO
391	!$OMP END PARALLEL
392
393	ENDIF
394
395	!
396	!-- Horizontal boundary conditions
397	CALL exchange_horiz( r, 1)
398
399	IF ( .NOT. bc_lr_cyc ) THEN
400	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
401	r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
402	ENDIF
403	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
404	r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
405	ENDIF
406	ENDIF
407
408	IF ( .NOT. bc_ns_cyc ) THEN
409	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
410	r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
411	ENDIF
412	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
413	r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
414	ENDIF
415	ENDIF
416
417	!
418	!-- Boundary conditions at bottom and top of the domain.outflow_l, outflow_n,
419	!-- These points are not handled by the above loop. Points may be within
420	!-- buildings, but that doesn't matter.
421	IF ( ibc_p_b == 1 ) THEN
422	r(nzb,:,: ) = r(nzb+1,:,:)
423	ELSE
424	r(nzb,:,: ) = 0.0_wp
425	ENDIF
426
427	IF ( ibc_p_t == 1 ) THEN
428	r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
429	ELSE
430	r(nzt_mg(l)+1,:,: ) = 0.0_wp
431	ENDIF
432
433	CALL cpu_log( log_point_s(53), 'resid', 'stop' )
434
435	END SUBROUTINE resid_fast
436
437
438	!------------------------------------------------------------------------------!
439	! Description:
440	! ------------
441	!> Interpolates the residual on the next coarser grid with "full weighting"
442	!> scheme
443	!------------------------------------------------------------------------------!
444	SUBROUTINE restrict_fast( f_mg, r )
445
446
447	USE control_parameters, &
448	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
449	inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l, &
450	nest_bound_n, nest_bound_r, nest_bound_s, outflow_l, &
451	outflow_n, outflow_r, outflow_s, topography
452
453	USE indices, &
454	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, &
455	wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, &
456	wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, &
457	nys_mg, nyn_mg, nzb, nzt_mg
458
459	IMPLICIT NONE
460
461	INTEGER(iwp) :: i !<
462	INTEGER(iwp) :: ic !<
463	INTEGER(iwp) :: j !<
464	INTEGER(iwp) :: jc !<
465	INTEGER(iwp) :: k !<
466	INTEGER(iwp) :: kc !<
467	INTEGER(iwp) :: l !<
468	INTEGER(iwp) :: km1 !<
469	INTEGER(iwp) :: kp1 !<
470
471	REAL(wp) :: rkjim !<
472	REAL(wp) :: rkjip !<
473	REAL(wp) :: rkjmi !<
474	REAL(wp) :: rkjmim !<
475	REAL(wp) :: rkjmip !<
476	REAL(wp) :: rkjpi !<
477	REAL(wp) :: rkjpim !<
478	REAL(wp) :: rkjpip !<
479	REAL(wp) :: rkmji !<
480	REAL(wp) :: rkmjim !<
481	REAL(wp) :: rkmjip !<
482	REAL(wp) :: rkmjmi !<
483	REAL(wp) :: rkmjmim !<
484	REAL(wp) :: rkmjmip !<
485	REAL(wp) :: rkmjpi !<
486	REAL(wp) :: rkmjpim !<
487	REAL(wp) :: rkmjpip !<
488
489	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
490	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
491	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
492
493	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1, &
494	nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
495	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !<
496
497	!
498	!-- Interpolate the residual
499	l = grid_level
500
501	CALL cpu_log( log_point_s(54), 'restrict', 'start' )
502
503	IF ( topography == 'flat' .OR. masking_method ) THEN
504	!
505	!-- No wall treatment
506	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
507	DO ic = nxl_mg(l), nxr_mg(l)
508	i = 2*ic
509	!$OMP DO SCHEDULE( STATIC )
510	DO jc = nys_mg(l), nyn_mg(l)
511	!
512	!-- Calculation for the first point along k
513	j = 2*jc
514	k = ind_even_odd+1
515	kp1 = k-ind_even_odd
516	kc = k-ind_even_odd
517	f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
518	8.0_wp * r(k,j,i) &
519	+ 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + &
520	r(k,j+1,i) + r(k,j-1,i) ) &
521	+ 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
522	r(k,j-1,i+1) + r(k,j+1,i+1) ) &
523	+ 16.0_wp * r(k,j,i) &
524	+ 4.0_wp * r(kp1,j,i) &
525	+ 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + &
526	r(kp1,j+1,i) + r(kp1,j-1,i) ) &
527	+ ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
528	r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
529	)
530	!
531	!-- Calculation for the other points along k
532	!DIR$ IVDEP
533	DO k = ind_even_odd+2, nzt_mg(l+1) ! Fine grid at this point
534	km1 = k-ind_even_odd-1
535	kp1 = k-ind_even_odd
536	kc = k-ind_even_odd ! Coarse grid index
537
538	f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
539	8.0_wp * r(k,j,i) &
540	+ 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + &
541	r(k,j+1,i) + r(k,j-1,i) ) &
542	+ 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
543	r(k,j-1,i+1) + r(k,j+1,i+1) ) &
544	+ 4.0_wp * r(km1,j,i) &
545	+ 2.0_wp * ( r(km1,j,i-1) + r(km1,j,i+1) + &
546	r(km1,j+1,i) + r(km1,j-1,i) ) &
547	+ ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
548	r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
549	+ 4.0_wp * r(kp1,j,i) &
550	+ 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + &
551	r(kp1,j+1,i) + r(kp1,j-1,i) ) &
552	+ ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
553	r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
554	)
555	ENDDO
556	ENDDO
557	!$OMP ENDDO nowait
558	ENDDO
559	!$OMP END PARALLEL
560
561	ELSE
562	!
563	!-- Choose flag array of the upper level
564	SELECT CASE ( l+1 )
565	CASE ( 1 )
566	flags => wall_flags_1
567	CASE ( 2 )
568	flags => wall_flags_2
569	CASE ( 3 )
570	flags => wall_flags_3
571	CASE ( 4 )
572	flags => wall_flags_4
573	CASE ( 5 )
574	flags => wall_flags_5
575	CASE ( 6 )
576	flags => wall_flags_6
577	CASE ( 7 )
578	flags => wall_flags_7
579	CASE ( 8 )
580	flags => wall_flags_8
581	CASE ( 9 )
582	flags => wall_flags_9
583	CASE ( 10 )
584	flags => wall_flags_10
585	END SELECT
586
587	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1,rkjim,rkjip,rkjpi, &
588	!$OMP rkjmim,rkjpim,rkjmip,rkjpip,rkmji,rkmjim, &
589	!$OMP rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim, &
590	!$OMP rkmjmip,rkmjpip)
591	!$OMP DO
592	DO ic = nxl_mg(l), nxr_mg(l)
593	i = 2*ic
594	DO jc = nys_mg(l), nyn_mg(l)
595	j = 2*jc
596	!DIR$ IVDEP
597	DO k = ind_even_odd+1, nzt_mg(l+1) !fine grid at this point
598
599	km1 = k-ind_even_odd-1
600	kp1 = k-ind_even_odd
601	kc = k-ind_even_odd !Coarse grid index
602	!
603	!-- Use implicit Neumann BCs if the respective gridpoint is
604	!-- inside the building
605	rkjim = MERGE(r(k,j,i),r(k,j,i-1),BTEST( flags(k,j,i-1), 6))
606	rkjip = MERGE(r(k,j,i),r(k,j,i+1),BTEST( flags(k,j,i+1), 6))
607	rkjpi = MERGE(r(k,j,i),r(k,j+1,i),BTEST( flags(k,j+1,i), 6))
608	rkjmi = MERGE(r(k,j,i),r(k,j-1,i),BTEST( flags(k,j-1,i), 6))
609
610	rkjmim = MERGE(r(k,j,i),r(k,j-1,i-1),BTEST( flags(k,j-1,i-1), 6))
611	rkjpim = MERGE(r(k,j,i),r(k,j+1,i-1),BTEST( flags(k,j+1,i-1), 6))
612	rkjmip = MERGE(r(k,j,i),r(k,j-1,i+1),BTEST( flags(k,j-1,i+1), 6))
613	rkjpip = MERGE(r(k,j,i),r(k,j+1,i+1),BTEST( flags(k,j+1,i+1), 6))
614
615	rkmji = MERGE(r(k,j,i),r(km1,j,i) ,BTEST( flags(km1,j,i) , 6))
616	rkmjim = MERGE(r(k,j,i),r(km1,j,i-1),BTEST( flags(km1,j,i-1), 6))
617	rkmjip = MERGE(r(k,j,i),r(km1,j,i+1),BTEST( flags(km1,j,i+1), 6))
618	rkmjpi = MERGE(r(k,j,i),r(km1,j+1,i),BTEST( flags(km1,j+1,i), 6))
619	rkmjmi = MERGE(r(k,j,i),r(km1,j-1,i),BTEST( flags(km1,j-1,i), 6))
620
621	rkmjmim = MERGE(r(k,j,i),r(km1,j-1,i-1),BTEST( flags(km1,j-1,i-1), 6))
622	rkmjpim = MERGE(r(k,j,i),r(km1,j+1,i-1),BTEST( flags(km1,j+1,i-1), 6))
623	rkmjmip = MERGE(r(k,j,i),r(km1,j-1,i+1),BTEST( flags(km1,j-1,i+1), 6))
624	rkmjpip = MERGE(r(k,j,i),r(km1,j+1,i+1),BTEST( flags(km1,j+1,i+1), 6))
625
626	f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
627	8.0_wp * r(k,j,i) &
628	+ 4.0_wp * ( rkjim + rkjip + rkjpi + rkjmi ) &
629	+ 2.0_wp * ( rkjmim + rkjpim + rkjmip + rkjpip ) &
630	+ 4.0_wp * rkmji &
631	+ 2.0_wp * ( rkmjim + rkmjip + rkmjpi + rkmjmi ) &
632	+ ( rkmjmim + rkmjpim + rkmjmip + rkmjpip ) &
633	+ 4.0_wp * r(kp1,j,i) &
634	+ 2.0_wp * ( r(kp1,j,i-1) + r(kp1,j,i+1) + &
635	r(kp1,j+1,i) + r(kp1,j-1,i) ) &
636	+ ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
637	r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
638	)
639	ENDDO
640	ENDDO
641	ENDDO
642	!$OMP END PARALLEL
643
644	ENDIF
645
646	!
647	!-- Horizontal boundary conditions
648	CALL exchange_horiz( f_mg, 1)
649
650	IF ( .NOT. bc_lr_cyc ) THEN
651	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
652	f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
653	ENDIF
654	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
655	f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
656	ENDIF
657	ENDIF
658
659	IF ( .NOT. bc_ns_cyc ) THEN
660	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
661	f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
662	ENDIF
663	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
664	f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
665	ENDIF
666	ENDIF
667
668	!
669	!-- Boundary conditions at bottom and top of the domain.
670	!-- These points are not handled by the above loop. Points may be within
671	!-- buildings, but that doesn't matter.
672	IF ( ibc_p_b == 1 ) THEN
673	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
674	ELSE
675	f_mg(nzb,:,: ) = 0.0_wp
676	ENDIF
677
678	IF ( ibc_p_t == 1 ) THEN
679	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
680	ELSE
681	f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
682	ENDIF
683
684	CALL cpu_log( log_point_s(54), 'restrict', 'stop' )
685	!
686	!-- Why do we need a sorting here?
687	CALL sort_k_to_even_odd_blocks( f_mg , l)
688
689	END SUBROUTINE restrict_fast
690
691
692	!------------------------------------------------------------------------------!
693	! Description:
694	! ------------
695	!> Interpolates the correction of the perturbation pressure
696	!> to the next finer grid.
697	!------------------------------------------------------------------------------!
698	SUBROUTINE prolong_fast( p, temp )
699
700
701	USE control_parameters, &
702	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
703	inflow_n, inflow_r, inflow_s, nest_bound_l, nest_bound_n, &
704	nest_bound_r, nest_bound_s, outflow_l, outflow_n, &
705	outflow_r, outflow_s
706
707	USE indices, &
708	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
709
710	IMPLICIT NONE
711
712	INTEGER(iwp) :: i !<
713	INTEGER(iwp) :: j !<
714	INTEGER(iwp) :: k !<
715	INTEGER(iwp) :: l !<
716	INTEGER(iwp) :: kp1 !<
717	INTEGER(iwp) :: ke !< Index for prolog even
718	INTEGER(iwp) :: ko !< Index for prolog odd
719
720	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
721	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
722	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !<
723
724	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
725	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
726	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !<
727
728
729	CALL cpu_log( log_point_s(55), 'prolong', 'start' )
730
731	!
732	!-- First, store elements of the coarser grid on the next finer grid
733	l = grid_level
734	ind_even_odd = even_odd_level(grid_level-1)
735
736	!$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko)
737	!$OMP DO
738	DO i = nxl_mg(l-1), nxr_mg(l-1)
739	DO j = nys_mg(l-1), nyn_mg(l-1)
740
741	!DIR$ IVDEP
742	DO k = ind_even_odd+1, nzt_mg(l-1)
743	kp1 = k - ind_even_odd
744	ke = 2 * ( k-ind_even_odd - 1 ) + 1
745	ko = 2 * k - 1
746	!
747	!-- Points of the coarse grid are directly stored on the next finer
748	!-- grid
749	temp(ko,2j,2i) = p(k,j,i)
750	!
751	!-- Points between two coarse-grid points
752	temp(ko,2j,2i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
753	temp(ko,2j+1,2i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
754	temp(ke,2j,2i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
755	!
756	!-- Points in the center of the planes stretched by four points
757	!-- of the coarse grid cube
758	temp(ko,2j+1,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
759	p(k,j+1,i) + p(k,j+1,i+1) )
760	temp(ke,2j,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
761	p(kp1,j,i) + p(kp1,j,i+1) )
762	temp(ke,2j+1,2i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + &
763	p(kp1,j,i) + p(kp1,j+1,i) )
764	!
765	!-- Points in the middle of coarse grid cube
766	temp(ke,2j+1,2i+1) = 0.125_wp * &
767	( p(k,j,i) + p(k,j,i+1) + &
768	p(k,j+1,i) + p(k,j+1,i+1) + &
769	p(kp1,j,i) + p(kp1,j,i+1) + &
770	p(kp1,j+1,i) + p(kp1,j+1,i+1) )
771	ENDDO
772
773	!DIR$ IVDEP
774	DO k = nzb+1, ind_even_odd
775	kp1 = k + ind_even_odd + 1
776	ke = 2 * k
777	ko = 2 * ( k + ind_even_odd )
778	!
779	!-- Points of the coarse grid are directly stored on the next finer
780	!-- grid
781	temp(ko,2j,2i) = p(k,j,i)
782	!
783	!-- Points between two coarse-grid points
784	temp(ko,2j,2i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
785	temp(ko,2j+1,2i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
786	temp(ke,2j,2i) = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
787	!
788	!-- Points in the center of the planes stretched by four points
789	!-- of the coarse grid cube
790	temp(ko,2j+1,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
791	p(k,j+1,i) + p(k,j+1,i+1) )
792	temp(ke,2j,2i+1) = 0.25_wp * ( p(k,j,i) + p(k,j,i+1) + &
793	p(kp1,j,i) + p(kp1,j,i+1) )
794	temp(ke,2j+1,2i) = 0.25_wp * ( p(k,j,i) + p(k,j+1,i) + &
795	p(kp1,j,i) + p(kp1,j+1,i) )
796	!
797	!-- Points in the middle of coarse grid cube
798	temp(ke,2j+1,2i+1) = 0.125_wp * &
799	( p(k,j,i) + p(k,j,i+1) + &
800	p(k,j+1,i) + p(k,j+1,i+1) + &
801	p(kp1,j,i) + p(kp1,j,i+1) + &
802	p(kp1,j+1,i) + p(kp1,j+1,i+1) )
803	ENDDO
804
805	ENDDO
806	ENDDO
807	!$OMP END PARALLEL
808
809	ind_even_odd = even_odd_level(grid_level)
810	!
811	!-- Horizontal boundary conditions
812	CALL exchange_horiz( temp, 1)
813
814	IF ( .NOT. bc_lr_cyc ) THEN
815	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
816	temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
817	ENDIF
818	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
819	temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
820	ENDIF
821	ENDIF
822
823	IF ( .NOT. bc_ns_cyc ) THEN
824	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
825	temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
826	ENDIF
827	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
828	temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
829	ENDIF
830	ENDIF
831
832	!
833	!-- Bottom and top boundary conditions
834	IF ( ibc_p_b == 1 ) THEN
835	temp(nzb,:,: ) = temp(ind_even_odd+1,:,:)
836	ELSE
837	temp(nzb,:,: ) = 0.0_wp
838	ENDIF
839
840	IF ( ibc_p_t == 1 ) THEN
841	temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:)
842	ELSE
843	temp(nzt_mg(l)+1,:,: ) = 0.0_wp
844	ENDIF
845
846	CALL cpu_log( log_point_s(55), 'prolong', 'stop' )
847
848	END SUBROUTINE prolong_fast
849
850
851	!------------------------------------------------------------------------------!
852	! Description:
853	! ------------
854	!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
855	!> 3D-Red-Black decomposition (GS-RB) is used.
856	!------------------------------------------------------------------------------!
857	SUBROUTINE redblack_fast( f_mg, p_mg )
858
859
860	USE arrays_3d, &
861	ONLY: f1_mg, f2_mg, f3_mg
862
863	USE control_parameters, &
864	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
865	inflow_n, inflow_r, inflow_s, masking_method, nest_bound_l, &
866	nest_bound_n, nest_bound_r, nest_bound_s, ngsrb, &
867	outflow_l, outflow_n, outflow_r, outflow_s, topography
868
869	USE grid_variables, &
870	ONLY: ddx2_mg, ddy2_mg
871
872	USE indices, &
873	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, &
874	wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, &
875	wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, &
876	nys_mg, nyn_mg, nzb, nzt_mg
877
878	IMPLICIT NONE
879
880	INTEGER(iwp) :: color !<
881	INTEGER(iwp) :: i !<
882	INTEGER(iwp) :: ic !<
883	INTEGER(iwp) :: j !<
884	INTEGER(iwp) :: jc !<
885	INTEGER(iwp) :: jj !<
886	INTEGER(iwp) :: k !<
887	INTEGER(iwp) :: l !<
888	INTEGER(iwp) :: n !<
889	INTEGER(iwp) :: km1 !<
890	INTEGER(iwp) :: kp1 !<
891
892	LOGICAL :: unroll !<
893
894	REAL(wp) :: wall_left !<
895	REAL(wp) :: wall_north !<
896	REAL(wp) :: wall_right !<
897	REAL(wp) :: wall_south !<
898	REAL(wp) :: wall_total !<
899	REAL(wp) :: wall_top !<
900
901	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
902	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
903	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
904	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
905	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
906	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
907
908	l = grid_level
909
910	!
911	!-- Choose flag array of this level
912	IF ( topography /= 'flat' .AND. .NOT. masking_method ) THEN
913
914	SELECT CASE ( l )
915	CASE ( 1 )
916	flags => wall_flags_1
917	CASE ( 2 )
918	flags => wall_flags_2
919	CASE ( 3 )
920	flags => wall_flags_3
921	CASE ( 4 )
922	flags => wall_flags_4
923	CASE ( 5 )
924	flags => wall_flags_5
925	CASE ( 6 )
926	flags => wall_flags_6
927	CASE ( 7 )
928	flags => wall_flags_7
929	CASE ( 8 )
930	flags => wall_flags_8
931	CASE ( 9 )
932	flags => wall_flags_9
933	CASE ( 10 )
934	flags => wall_flags_10
935	END SELECT
936
937	ENDIF
938
939	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
940	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
941
942	DO n = 1, ngsrb
943
944	DO color = 1, 2
945	!
946	!-- No wall treatment in case of flat topography
947	IF ( topography == 'flat' .OR. masking_method ) THEN
948
949	IF ( .NOT. unroll ) THEN
950
951	CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
952	!
953	!-- Without unrolling of loops, no cache optimization
954	!$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
955	!$OMP DO
956	DO i = nxl_mg(l), nxr_mg(l), 2
957	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
958	!DIR$ IVDEP
959	DO k = ind_even_odd+1, nzt_mg(l)
960	km1 = k-ind_even_odd-1
961	kp1 = k-ind_even_odd
962	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
963	ddx2_mg(l) * &
964	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
965	+ ddy2_mg(l) * &
966	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
967	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
968	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
969	- f_mg(k,j,i) )
970	ENDDO
971	ENDDO
972	ENDDO
973
974	!$OMP DO
975	DO i = nxl_mg(l)+1, nxr_mg(l), 2
976	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
977	!DIR$ IVDEP
978	DO k = ind_even_odd+1, nzt_mg(l)
979	km1 = k-ind_even_odd-1
980	kp1 = k-ind_even_odd
981	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
982	ddx2_mg(l) * &
983	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
984	+ ddy2_mg(l) * &
985	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
986	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
987	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
988	- f_mg(k,j,i) )
989	ENDDO
990	ENDDO
991	ENDDO
992
993	!$OMP DO
994	DO i = nxl_mg(l), nxr_mg(l), 2
995	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
996	!DIR$ IVDEP
997	DO k = nzb+1, ind_even_odd
998	km1 = k+ind_even_odd
999	kp1 = k+ind_even_odd+1
1000	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1001	ddx2_mg(l) * &
1002	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1003	+ ddy2_mg(l) * &
1004	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1005	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1006	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1007	- f_mg(k,j,i) )
1008	ENDDO
1009	ENDDO
1010	ENDDO
1011
1012	!$OMP DO
1013	DO i = nxl_mg(l)+1, nxr_mg(l), 2
1014	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1015	!DIR$ IVDEP
1016	DO k = nzb+1, ind_even_odd
1017	km1 = k+ind_even_odd
1018	kp1 = k+ind_even_odd+1
1019	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1020	ddx2_mg(l) * &
1021	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1022	+ ddy2_mg(l) * &
1023	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1024	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1025	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1026	- f_mg(k,j,i) )
1027	ENDDO
1028	ENDDO
1029	ENDDO
1030	!$OMP END PARALLEL
1031
1032	CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
1033
1034	ELSE
1035	!
1036	!-- Loop unrolling along y, only one i loop for better cache use
1037	CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
1038
1039	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
1040	!$OMP DO
1041	DO ic = nxl_mg(l), nxr_mg(l), 2
1042	DO jc = nys_mg(l), nyn_mg(l), 4
1043	i = ic
1044	jj = jc+2-color
1045	!DIR$ IVDEP
1046	DO k = ind_even_odd+1, nzt_mg(l)
1047	km1 = k-ind_even_odd-1
1048	kp1 = k-ind_even_odd
1049	j = jj
1050	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1051	ddx2_mg(l) * &
1052	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1053	+ ddy2_mg(l) * &
1054	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1055	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1056	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1057	- f_mg(k,j,i) )
1058	j = jj+2
1059	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1060	ddx2_mg(l) * &
1061	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1062	+ ddy2_mg(l) * &
1063	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1064	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1065	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1066	- f_mg(k,j,i) )
1067	ENDDO
1068
1069	i = ic+1
1070	jj = jc+color-1
1071	!DIR$ IVDEP
1072	DO k = ind_even_odd+1, nzt_mg(l)
1073	km1 = k-ind_even_odd-1
1074	kp1 = k-ind_even_odd
1075	j = jj
1076	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1077	ddx2_mg(l) * &
1078	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1079	+ ddy2_mg(l) * &
1080	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1081	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1082	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1083	- f_mg(k,j,i) )
1084	j = jj+2
1085	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1086	ddx2_mg(l) * &
1087	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1088	+ ddy2_mg(l) * &
1089	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1090	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1091	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1092	- f_mg(k,j,i) )
1093	ENDDO
1094
1095	i = ic
1096	jj = jc+color-1
1097	!DIR$ IVDEP
1098	DO k = nzb+1, ind_even_odd
1099	km1 = k+ind_even_odd
1100	kp1 = k+ind_even_odd+1
1101	j = jj
1102	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1103	ddx2_mg(l) * &
1104	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1105	+ ddy2_mg(l) * &
1106	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1107	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1108	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1109	- f_mg(k,j,i) )
1110	j = jj+2
1111	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1112	ddx2_mg(l) * &
1113	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1114	+ ddy2_mg(l) * &
1115	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1116	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1117	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1118	- f_mg(k,j,i) )
1119	ENDDO
1120
1121	i = ic+1
1122	jj = jc+2-color
1123	!DIR$ IVDEP
1124	DO k = nzb+1, ind_even_odd
1125	km1 = k+ind_even_odd
1126	kp1 = k+ind_even_odd+1
1127	j = jj
1128	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1129	ddx2_mg(l) * &
1130	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1131	+ ddy2_mg(l) * &
1132	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1133	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1134	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1135	- f_mg(k,j,i) )
1136	j = jj+2
1137	p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * ( &
1138	ddx2_mg(l) * &
1139	( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
1140	+ ddy2_mg(l) * &
1141	( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
1142	+ f2_mg_b(k,l) * p_mg(kp1,j,i) &
1143	+ f3_mg_b(k,l) * p_mg(km1,j,i) &
1144	- f_mg(k,j,i) )
1145	ENDDO
1146
1147	ENDDO
1148	ENDDO
1149	!$OMP END PARALLEL
1150
1151	CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
1152
1153	ENDIF
1154
1155	ELSE
1156
1157	IF ( .NOT. unroll ) THEN
1158
1159	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
1160	!
1161	!-- Without unrolling of loops, no cache optimization /
1162	!-- vectorization. Therefore, the non-vectorized IBITS function
1163	!-- is still used.
1164	!$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
1165	!$OMP DO
1166	DO i = nxl_mg(l), nxr_mg(l), 2
1167	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1168	!DIR$ IVDEP
1169	DO k = ind_even_odd+1, nzt_mg(l)
1170	km1 = k-ind_even_odd-1
1171	kp1 = k-ind_even_odd
1172
1173	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
1174	ddx2_mg(l) * &
1175	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1176	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1177	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1178	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1179	+ ddy2_mg(l) * &
1180	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1181	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1182	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1183	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1184	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1185	+ f3_mg(k,l) * &
1186	( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1187	( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
1188	- f_mg(k,j,i) )
1189	ENDDO
1190	ENDDO
1191	ENDDO
1192
1193	!$OMP DO
1194	DO i = nxl_mg(l)+1, nxr_mg(l), 2
1195	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
1196	!DIR$ IVDEP
1197	DO k = ind_even_odd+1, nzt_mg(l)
1198	km1 = k-ind_even_odd-1
1199	kp1 = k-ind_even_odd
1200
1201	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
1202	ddx2_mg(l) * &
1203	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1204	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1205	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1206	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1207	+ ddy2_mg(l) * &
1208	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1209	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1210	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1211	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1212	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1213	+ f3_mg(k,l) * &
1214	( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1215	( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
1216	- f_mg(k,j,i) )
1217	ENDDO
1218	ENDDO
1219	ENDDO
1220
1221	!$OMP DO
1222	DO i = nxl_mg(l), nxr_mg(l), 2
1223	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
1224	!DIR$ IVDEP
1225	DO k = nzb+1, ind_even_odd
1226	km1 = k+ind_even_odd
1227	kp1 = k+ind_even_odd+1
1228	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
1229	ddx2_mg(l) * &
1230	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1231	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1232	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1233	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1234	+ ddy2_mg(l) * &
1235	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1236	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1237	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1238	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1239	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1240	+ f3_mg(k,l) * &
1241	( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1242	( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
1243	- f_mg(k,j,i) )
1244	ENDDO
1245	ENDDO
1246	ENDDO
1247
1248	!$OMP DO
1249	DO i = nxl_mg(l)+1, nxr_mg(l), 2
1250	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1251	!DIR$ IVDEP
1252	DO k = nzb+1, ind_even_odd
1253	km1 = k+ind_even_odd
1254	kp1 = k+ind_even_odd+1
1255
1256	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
1257	ddx2_mg(l) * &
1258	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1259	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1260	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1261	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1262	+ ddy2_mg(l) * &
1263	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1264	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1265	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1266	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1267	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1268	+ f3_mg(k,l) * &
1269	( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1270	( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
1271	- f_mg(k,j,i) )
1272	ENDDO
1273	ENDDO
1274	ENDDO
1275	!$OMP END PARALLEL
1276
1277	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
1278
1279	ELSE
1280	!
1281	!-- Loop unrolling along y, only one i loop for better cache use
1282	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
1283
1284	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
1285	!$OMP DO
1286	DO ic = nxl_mg(l), nxr_mg(l), 2
1287	DO jc = nys_mg(l), nyn_mg(l), 4
1288	i = ic
1289	jj = jc+2-color
1290	!DIR$ IVDEP
1291	DO k = ind_even_odd+1, nzt_mg(l)
1292
1293	km1 = k-ind_even_odd-1
1294	kp1 = k-ind_even_odd
1295
1296	j = jj
1297	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1298	ddx2_mg(l) * &
1299	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1300	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1301	+ ddy2_mg(l) * &
1302	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1303	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1304	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1305	+ f3_mg(k,l) * &
1306	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1307	- f_mg(k,j,i) )
1308	j = jj+2
1309	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1310	ddx2_mg(l) * &
1311	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1312	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1313	+ ddy2_mg(l) * &
1314	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1315	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1316	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1317	+ f3_mg(k,l) * &
1318	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1319	- f_mg(k,j,i) )
1320
1321	ENDDO
1322
1323	i = ic+1
1324	jj = jc+color-1
1325	!DIR$ IVDEP
1326	DO k = ind_even_odd+1, nzt_mg(l)
1327
1328	km1 = k-ind_even_odd-1
1329	kp1 = k-ind_even_odd
1330
1331	j =jj
1332	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1333	ddx2_mg(l) * &
1334	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1335	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1336	+ ddy2_mg(l) * &
1337	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1338	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1339	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1340	+ f3_mg(k,l) * &
1341	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1342	- f_mg(k,j,i) )
1343
1344	j = jj+2
1345	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1346	ddx2_mg(l) * &
1347	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1348	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1349	+ ddy2_mg(l) * &
1350	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1351	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1352	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1353	+ f3_mg(k,l) * &
1354	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1355	- f_mg(k,j,i) )
1356
1357	ENDDO
1358
1359	i = ic
1360	jj = jc+color-1
1361	!DIR$ IVDEP
1362	DO k = nzb+1, ind_even_odd
1363
1364	km1 = k+ind_even_odd
1365	kp1 = k+ind_even_odd+1
1366
1367	j =jj
1368	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1369	ddx2_mg(l) * &
1370	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1371	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1372	+ ddy2_mg(l) * &
1373	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1374	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1375	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1376	+ f3_mg(k,l) * &
1377	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1378	- f_mg(k,j,i) )
1379
1380	j = jj+2
1381	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1382	ddx2_mg(l) * &
1383	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1384	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1385	+ ddy2_mg(l) * &
1386	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1387	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1388	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1389	+ f3_mg(k,l) * &
1390	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1391	- f_mg(k,j,i) )
1392
1393	ENDDO
1394
1395	i = ic+1
1396	jj = jc+2-color
1397	!DIR$ IVDEP
1398	DO k = nzb+1, ind_even_odd
1399
1400	km1 = k+ind_even_odd
1401	kp1 = k+ind_even_odd+1
1402
1403	j =jj
1404	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1405	ddx2_mg(l) * &
1406	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1407	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1408	+ ddy2_mg(l) * &
1409	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1410	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1411	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1412	+ f3_mg(k,l) * &
1413	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1414	- f_mg(k,j,i) )
1415
1416	j = jj+2
1417	p_mg(k,j,i) = 1.0 / f1_mg(k,l) * ( &
1418	ddx2_mg(l) * &
1419	( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5)) &
1420	+ MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1421	+ ddy2_mg(l) * &
1422	( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3)) &
1423	+ MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1424	+ f2_mg(k,l) * p_mg(kp1,j,i) &
1425	+ f3_mg(k,l) * &
1426	MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0)) &
1427	- f_mg(k,j,i) )
1428
1429	ENDDO
1430
1431	ENDDO
1432	ENDDO
1433	!$OMP END PARALLEL
1434
1435	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
1436
1437	ENDIF
1438
1439	ENDIF
1440
1441	!
1442	!-- Horizontal boundary conditions
1443	CALL special_exchange_horiz( p_mg, color )
1444
1445	IF ( .NOT. bc_lr_cyc ) THEN
1446	IF ( inflow_l .OR. outflow_l .OR. nest_bound_l ) THEN
1447	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
1448	ENDIF
1449	IF ( inflow_r .OR. outflow_r .OR. nest_bound_r ) THEN
1450	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
1451	ENDIF
1452	ENDIF
1453
1454	IF ( .NOT. bc_ns_cyc ) THEN
1455	IF ( inflow_n .OR. outflow_n .OR. nest_bound_n ) THEN
1456	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
1457	ENDIF
1458	IF ( inflow_s .OR. outflow_s .OR. nest_bound_s ) THEN
1459	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
1460	ENDIF
1461	ENDIF
1462
1463	!
1464	!-- Bottom and top boundary conditions
1465	IF ( ibc_p_b == 1 ) THEN
1466	p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:)
1467	ELSE
1468	p_mg(nzb,:,: ) = 0.0_wp
1469	ENDIF
1470
1471	IF ( ibc_p_t == 1 ) THEN
1472	p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:)
1473	ELSE
1474	p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
1475	ENDIF
1476
1477	ENDDO
1478	ENDDO
1479	!
1480	!-- Set pressure within topography and at the topography surfaces
1481	IF ( topography /= 'flat' .AND. .NOT. masking_method ) THEN
1482
1483	!$OMP PARALLEL PRIVATE (i,j,k,kp1,wall_left,wall_north,wall_right, &
1484	!$OMP wall_south,wall_top,wall_total)
1485	!$OMP DO
1486	DO i = nxl_mg(l), nxr_mg(l)
1487	DO j = nys_mg(l), nyn_mg(l)
1488	!DIR$ IVDEP
1489	DO k = ind_even_odd+1, nzt_mg(l)
1490	km1 = k-ind_even_odd-1
1491	kp1 = k-ind_even_odd
1492	!
1493	!-- First, set pressure inside topography to zero
1494	p_mg(k,j,i) = p_mg(k,j,i) * &
1495	( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1496	!
1497	!-- Second, determine if the gridpoint inside topography is
1498	!-- adjacent to a wall and set its value to a value given by the
1499	!-- average of those values obtained from Neumann boundary
1500	!-- condition
1501	wall_left = IBITS( flags(k,j,i-1), 5, 1 )
1502	wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1503	wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1504	wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1505	wall_top = IBITS( flags(kp1,j,i), 0, 1 )
1506	wall_total = wall_left + wall_right + wall_south + &
1507	wall_north + wall_top
1508
1509	IF ( wall_total > 0.0_wp ) THEN
1510	p_mg(k,j,i) = 1.0_wp / wall_total * &
1511	( wall_left * p_mg(k,j,i-1) + &
1512	wall_right * p_mg(k,j,i+1) + &
1513	wall_south * p_mg(k,j-1,i) + &
1514	wall_north * p_mg(k,j+1,i) + &
1515	wall_top * p_mg(kp1,j,i) )
1516	ENDIF
1517	ENDDO
1518
1519	!DIR$ IVDEP
1520	DO k = nzb+1, ind_even_odd
1521	km1 = k+ind_even_odd
1522	kp1 = k+ind_even_odd+1
1523	!
1524	!-- First, set pressure inside topography to zero
1525	p_mg(k,j,i) = p_mg(k,j,i) * &
1526	( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1527	!
1528	!-- Second, determine if the gridpoint inside topography is
1529	!-- adjacent to a wall and set its value to a value given by the
1530	!-- average of those values obtained from Neumann boundary
1531	!-- condition
1532	wall_left = IBITS( flags(k,j,i-1), 5, 1 )
1533	wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1534	wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1535	wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1536	wall_top = IBITS( flags(kp1,j,i), 0, 1 )
1537	wall_total = wall_left + wall_right + wall_south + &
1538	wall_north + wall_top
1539
1540	IF ( wall_total > 0.0_wp ) THEN
1541	p_mg(k,j,i) = 1.0_wp / wall_total * &
1542	( wall_left * p_mg(k,j,i-1) + &
1543	wall_right * p_mg(k,j,i+1) + &
1544	wall_south * p_mg(k,j-1,i) + &
1545	wall_north * p_mg(k,j+1,i) + &
1546	wall_top * p_mg(kp1,j,i) )
1547	ENDIF
1548	ENDDO
1549	ENDDO
1550	ENDDO
1551	!$OMP END PARALLEL
1552
1553	!
1554	!-- One more time horizontal boundary conditions
1555	CALL exchange_horiz( p_mg, 1)
1556
1557	ENDIF
1558
1559	END SUBROUTINE redblack_fast
1560
1561
1562	!------------------------------------------------------------------------------!
1563	! Description:
1564	! ------------
1565	!> Sort k-Dimension from sequential into blocks of even and odd.
1566	!> This is required to vectorize the red-black subroutine.
1567	!> Version for 3D-REAL arrays
1568	!------------------------------------------------------------------------------!
1569	SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel )
1570
1571
1572	USE control_parameters, &
1573	ONLY: grid_level
1574
1575	USE indices, &
1576	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1577
1578	IMPLICIT NONE
1579
1580	INTEGER(iwp), INTENT(IN) :: glevel
1581
1582	REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1, &
1583	nys_mg(glevel)-1:nyn_mg(glevel)+1, &
1584	nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: p_mg !<
1585	!
1586	!-- Local variables
1587	INTEGER(iwp) :: i !<
1588	INTEGER(iwp) :: j !<
1589	INTEGER(iwp) :: k !<
1590	INTEGER(iwp) :: l !<
1591	INTEGER(iwp) :: ind !<
1592	REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp !<
1593
1594
1595	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
1596
1597	l = glevel
1598	ind_even_odd = even_odd_level(l)
1599
1600	!$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
1601	!$OMP DO
1602	DO i = nxl_mg(l)-1, nxr_mg(l)+1
1603	DO j = nys_mg(l)-1, nyn_mg(l)+1
1604
1605	!
1606	!-- Sort the data with even k index
1607	ind = nzb-1
1608	DO k = nzb, nzt_mg(l), 2
1609	ind = ind + 1
1610	tmp(ind) = p_mg(k,j,i)
1611	ENDDO
1612	!
1613	!-- Sort the data with odd k index
1614	DO k = nzb+1, nzt_mg(l)+1, 2
1615	ind = ind + 1
1616	tmp(ind) = p_mg(k,j,i)
1617	ENDDO
1618
1619	p_mg(:,j,i) = tmp
1620
1621	ENDDO
1622	ENDDO
1623	!$OMP END PARALLEL
1624
1625	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
1626
1627	END SUBROUTINE sort_k_to_even_odd_blocks
1628
1629
1630	!------------------------------------------------------------------------------!
1631	! Description:
1632	! ------------
1633	!> Sort k-Dimension from sequential into blocks of even and odd.
1634	!> This is required to vectorize the red-black subroutine.
1635	!> Version for 1D-REAL arrays
1636	!------------------------------------------------------------------------------!
1637	SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel )
1638
1639
1640	USE indices, &
1641	ONLY: nzb, nzt_mg
1642
1643	IMPLICIT NONE
1644
1645	INTEGER(iwp), INTENT(IN) :: glevel
1646
1647	REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) :: f_mg
1648	REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) :: f_mg_b
1649
1650	!
1651	!-- Local variables
1652	INTEGER(iwp) :: ind !<
1653	INTEGER(iwp) :: k !<
1654
1655
1656	ind = nzb - 1
1657	!
1658	!-- Sort the data with even k index
1659	DO k = nzb, nzt_mg(glevel), 2
1660	ind = ind + 1
1661	IF ( k >= nzb+1 .AND. k <= nzt_mg(glevel) ) THEN
1662	f_mg_b(ind) = f_mg(k)
1663	ENDIF
1664	ENDDO
1665	!
1666	!-- Sort the data with odd k index
1667	DO k = nzb+1, nzt_mg(glevel)+1, 2
1668	ind = ind + 1
1669	IF( k >= nzb+1 .AND. k <= nzt_mg(glevel) ) THEN
1670	f_mg_b(ind) = f_mg(k)
1671	ENDIF
1672	ENDDO
1673
1674	END SUBROUTINE sort_k_to_even_odd_blocks_1d
1675
1676
1677	!------------------------------------------------------------------------------!
1678	! Description:
1679	! ------------
1680	!> Sort k-Dimension from sequential into blocks of even and odd.
1681	!> This is required to vectorize the red-black subroutine.
1682	!> Version for 2D-INTEGER arrays
1683	!------------------------------------------------------------------------------!
1684	SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
1685
1686
1687	USE control_parameters, &
1688	ONLY: grid_level
1689
1690	USE indices, &
1691	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1692
1693	IMPLICIT NONE
1694
1695	INTEGER(iwp), INTENT(IN) :: glevel
1696
1697	INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1, &
1698	nys_mg(glevel)-1:nyn_mg(glevel)+1, &
1699	nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: i_mg !<
1700	!
1701	!-- Local variables
1702	INTEGER(iwp) :: i !<
1703	INTEGER(iwp) :: j !<
1704	INTEGER(iwp) :: k !<
1705	INTEGER(iwp) :: l !<
1706	INTEGER(iwp) :: ind !<
1707	INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) :: tmp
1708
1709
1710	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
1711
1712	l = glevel
1713	ind_even_odd = even_odd_level(l)
1714
1715	DO i = nxl_mg(l)-1, nxr_mg(l)+1
1716	DO j = nys_mg(l)-1, nyn_mg(l)+1
1717
1718	!
1719	!-- Sort the data with even k index
1720	ind = nzb-1
1721	DO k = nzb, nzt_mg(l), 2
1722	ind = ind + 1
1723	tmp(ind) = i_mg(k,j,i)
1724	ENDDO
1725	!
1726	!++ ATTENTION: Check reason for this error. Remove it or replace WRITE
1727	!++ by PALM message
1728	#if defined ( __parallel )
1729	IF ( ind /= ind_even_odd ) THEN
1730	WRITE (0,*) 'ERROR ==> illegal ind_even_odd ',ind,ind_even_odd,l
1731	CALL MPI_ABORT(MPI_COMM_WORLD,i,j)
1732	ENDIF
1733	#endif
1734	!
1735	!-- Sort the data with odd k index
1736	DO k = nzb+1, nzt_mg(l)+1, 2
1737	ind = ind + 1
1738	tmp(ind) = i_mg(k,j,i)
1739	ENDDO
1740
1741	i_mg(:,j,i) = tmp
1742
1743	ENDDO
1744	ENDDO
1745
1746	CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
1747
1748	END SUBROUTINE sort_k_to_even_odd_blocks_int
1749
1750
1751	!------------------------------------------------------------------------------!
1752	! Description:
1753	! ------------
1754	!> Sort k-dimension from blocks of even and odd into sequential
1755	!------------------------------------------------------------------------------!
1756	SUBROUTINE sort_k_to_sequential( p_mg )
1757
1758
1759	USE control_parameters, &
1760	ONLY: grid_level
1761
1762	USE indices, &
1763	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1764
1765	IMPLICIT NONE
1766
1767	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1768	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1769	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
1770	!
1771	!-- Local variables
1772	INTEGER(iwp) :: i !<
1773	INTEGER(iwp) :: j !<
1774	INTEGER(iwp) :: k !<
1775	INTEGER(iwp) :: l !<
1776	INTEGER(iwp) :: ind !<
1777
1778	REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) :: tmp
1779
1780
1781	l = grid_level
1782
1783	!$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
1784	!$OMP DO
1785	DO i = nxl_mg(l)-1, nxr_mg(l)+1
1786	DO j = nys_mg(l)-1, nyn_mg(l)+1
1787
1788	ind = nzb - 1
1789	tmp = p_mg(:,j,i)
1790	DO k = nzb, nzt_mg(l), 2
1791	ind = ind + 1
1792	p_mg(k,j,i) = tmp(ind)
1793	ENDDO
1794
1795	DO k = nzb+1, nzt_mg(l)+1, 2
1796	ind = ind + 1
1797	p_mg(k,j,i) = tmp(ind)
1798	ENDDO
1799	ENDDO
1800	ENDDO
1801	!$OMP END PARALLEL
1802
1803	END SUBROUTINE sort_k_to_sequential
1804
1805
1806	!------------------------------------------------------------------------------!
1807	! Description:
1808	! ------------
1809	!> Gather subdomain data from all PEs.
1810	!------------------------------------------------------------------------------!
1811	SUBROUTINE mg_gather_fast( f2, f2_sub )
1812
1813	USE control_parameters, &
1814	ONLY: grid_level
1815
1816	USE cpulog, &
1817	ONLY: cpu_log, log_point_s
1818
1819	USE indices, &
1820	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1821
1822	IMPLICIT NONE
1823
1824	INTEGER(iwp) :: i !<
1825	INTEGER(iwp) :: il !<
1826	INTEGER(iwp) :: ir !<
1827	INTEGER(iwp) :: j !<
1828	INTEGER(iwp) :: jn !<
1829	INTEGER(iwp) :: js !<
1830	INTEGER(iwp) :: k !<
1831	INTEGER(iwp) :: nwords !<
1832
1833	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1834	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1835	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !<
1836	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1837	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1838	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !<
1839
1840	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
1841	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1842	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !<
1843
1844
1845	#if defined( __parallel )
1846	CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
1847
1848	f2_l = 0.0_wp
1849
1850	!
1851	!-- Store the local subdomain array on the total array
1852	js = mg_loc_ind(3,myid)
1853	IF ( south_border_pe ) js = js - 1
1854	jn = mg_loc_ind(4,myid)
1855	IF ( north_border_pe ) jn = jn + 1
1856	il = mg_loc_ind(1,myid)
1857	IF ( left_border_pe ) il = il - 1
1858	ir = mg_loc_ind(2,myid)
1859	IF ( right_border_pe ) ir = ir + 1
1860	DO i = il, ir
1861	DO j = js, jn
1862	DO k = nzb, nzt_mg(grid_level)+1
1863	f2_l(k,j,i) = f2_sub(k,j,i)
1864	ENDDO
1865	ENDDO
1866	ENDDO
1867
1868	!
1869	!-- Find out the number of array elements of the total array
1870	nwords = SIZE( f2 )
1871
1872	!
1873	!-- Gather subdomain data from all PEs
1874	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
1875	CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1876	f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1877	nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1878
1879	CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
1880	#endif
1881
1882	END SUBROUTINE mg_gather_fast
1883
1884
1885
1886	!------------------------------------------------------------------------------!
1887	! Description:
1888	! ------------
1889	!> @todo It might be possible to improve the speed of this routine by using
1890	!> non-blocking communication
1891	!------------------------------------------------------------------------------!
1892	SUBROUTINE mg_scatter_fast( p2, p2_sub )
1893
1894	USE control_parameters, &
1895	ONLY: grid_level
1896
1897	USE cpulog, &
1898	ONLY: cpu_log, log_point_s
1899
1900	USE indices, &
1901	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1902
1903	IMPLICIT NONE
1904
1905	INTEGER(iwp) :: nwords !<
1906
1907	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
1908	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
1909	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
1910
1911	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
1912	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1913	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !<
1914
1915	!
1916	!-- Find out the number of array elements of the subdomain array
1917	nwords = SIZE( p2_sub )
1918
1919	#if defined( __parallel )
1920	CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
1921
1922	p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1923	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1924
1925	CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
1926	#endif
1927
1928	END SUBROUTINE mg_scatter_fast
1929
1930
1931	!------------------------------------------------------------------------------!
1932	! Description:
1933	! ------------
1934	!> This is where the multigrid technique takes place. V- and W- Cycle are
1935	!> implemented and steered by the parameter "gamma". Parameter "nue" determines
1936	!> the convergence of the multigrid iterative solution. There are nue times
1937	!> RB-GS iterations. It should be set to "1" or "2", considering the time effort
1938	!> one would like to invest. Last choice shows a very good converging factor,
1939	!> but leads to an increase in computing time.
1940	!------------------------------------------------------------------------------!
1941	RECURSIVE SUBROUTINE next_mg_level_fast( f_mg, p_mg, p3, r )
1942
1943	USE control_parameters, &
1944	ONLY: bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir, &
1945	gamma_mg, grid_level, grid_level_count, ibc_p_b, ibc_p_t, &
1946	inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level, &
1947	mg_switch_to_pe0_level, mg_switch_to_pe0, nest_domain, &
1948	nest_bound_l, nest_bound_n, nest_bound_r, nest_bound_s, &
1949	ngsrb, outflow_l, outflow_n, outflow_r, outflow_s
1950
1951	USE indices, &
1952	ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, &
1953	nyn_mg, nzb, nzt, nzt_mg
1954
1955	IMPLICIT NONE
1956
1957	INTEGER(iwp) :: i !<
1958	INTEGER(iwp) :: j !<
1959	INTEGER(iwp) :: k !<
1960	INTEGER(iwp) :: nxl_mg_save !<
1961	INTEGER(iwp) :: nxr_mg_save !<
1962	INTEGER(iwp) :: nyn_mg_save !<
1963	INTEGER(iwp) :: nys_mg_save !<
1964	INTEGER(iwp) :: nzt_mg_save !<
1965
1966	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1967	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1968	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
1969	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1970	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1971	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
1972	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1973	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1974	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !<
1975	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
1976	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
1977	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
1978
1979	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
1980	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
1981	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !<
1982	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
1983	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
1984	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
1985
1986	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !<
1987	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !<
1988
1989	!
1990	!-- Restriction to the coarsest grid
1991	10 IF ( grid_level == 1 ) THEN
1992
1993	!
1994	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel
1995	!-- iterations in order to get a more accurate solution.
1996	ngsrb = 2 * ngsrb
1997
1998	ind_even_odd = even_odd_level(grid_level)
1999
2000	CALL redblack_fast( f_mg, p_mg )
2001
2002	ngsrb = ngsrb / 2
2003
2004
2005	ELSEIF ( grid_level /= 1 ) THEN
2006
2007	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
2008
2009	!
2010	!-- Solution on the actual grid level
2011	ind_even_odd = even_odd_level(grid_level)
2012
2013	CALL redblack_fast( f_mg, p_mg )
2014
2015	!
2016	!-- Determination of the actual residual
2017	CALL resid_fast( f_mg, p_mg, r )
2018
2019	!-- Restriction of the residual (finer grid values!) to the next coarser
2020	!-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have
2021	!-- to be set to the coarse grid values, because these variables are needed
2022	!-- for the exchange of ghost points in routine exchange_horiz
2023	grid_level = grid_level - 1
2024
2025	nxl = nxl_mg(grid_level)
2026	nys = nys_mg(grid_level)
2027	nxr = nxr_mg(grid_level)
2028	nyn = nyn_mg(grid_level)
2029	nzt = nzt_mg(grid_level)
2030
2031	IF ( grid_level == mg_switch_to_pe0_level ) THEN
2032
2033	!
2034	!-- From this level on, calculations are done on PE0 only.
2035	!-- First, carry out restriction on the subdomain.
2036	!-- Therefore, indices of the level have to be changed to subdomain
2037	!-- values in between (otherwise, the restrict routine would expect
2038	!-- the gathered array)
2039
2040	nxl_mg_save = nxl_mg(grid_level)
2041	nxr_mg_save = nxr_mg(grid_level)
2042	nys_mg_save = nys_mg(grid_level)
2043	nyn_mg_save = nyn_mg(grid_level)
2044	nzt_mg_save = nzt_mg(grid_level)
2045	nxl_mg(grid_level) = mg_loc_ind(1,myid)
2046	nxr_mg(grid_level) = mg_loc_ind(2,myid)
2047	nys_mg(grid_level) = mg_loc_ind(3,myid)
2048	nyn_mg(grid_level) = mg_loc_ind(4,myid)
2049	nzt_mg(grid_level) = mg_loc_ind(5,myid)
2050	nxl = mg_loc_ind(1,myid)
2051	nxr = mg_loc_ind(2,myid)
2052	nys = mg_loc_ind(3,myid)
2053	nyn = mg_loc_ind(4,myid)
2054	nzt = mg_loc_ind(5,myid)
2055
2056	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, &
2057	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
2058	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
2059
2060	CALL restrict_fast( f2_sub, r )
2061
2062	!
2063	!-- Restore the correct indices of this level
2064	nxl_mg(grid_level) = nxl_mg_save
2065	nxr_mg(grid_level) = nxr_mg_save
2066	nys_mg(grid_level) = nys_mg_save
2067	nyn_mg(grid_level) = nyn_mg_save
2068	nzt_mg(grid_level) = nzt_mg_save
2069	nxl = nxl_mg(grid_level)
2070	nxr = nxr_mg(grid_level)
2071	nys = nys_mg(grid_level)
2072	nyn = nyn_mg(grid_level)
2073	nzt = nzt_mg(grid_level)
2074	!
2075	!-- Gather all arrays from the subdomains on PE0
2076	CALL mg_gather_fast( f2, f2_sub )
2077
2078	!
2079	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange
2080	!-- has to be carried out from now on
2081	mg_switch_to_pe0 = .TRUE.
2082
2083	!
2084	!-- In case of non-cyclic lateral boundary conditions, both in- and
2085	!-- outflow conditions have to be used on all PEs after the switch,
2086	!-- because then they have the total domain.
2087	IF ( bc_lr_dirrad ) THEN
2088	inflow_l = .TRUE.
2089	inflow_r = .FALSE.
2090	outflow_l = .FALSE.
2091	outflow_r = .TRUE.
2092	ELSEIF ( bc_lr_raddir ) THEN
2093	inflow_l = .FALSE.
2094	inflow_r = .TRUE.
2095	outflow_l = .TRUE.
2096	outflow_r = .FALSE.
2097	ELSEIF ( nest_domain ) THEN
2098	nest_bound_l = .TRUE.
2099	nest_bound_r = .TRUE.
2100	ENDIF
2101
2102	IF ( bc_ns_dirrad ) THEN
2103	inflow_n = .TRUE.
2104	inflow_s = .FALSE.
2105	outflow_n = .FALSE.
2106	outflow_s = .TRUE.
2107	ELSEIF ( bc_ns_raddir ) THEN
2108	inflow_n = .FALSE.
2109	inflow_s = .TRUE.
2110	outflow_n = .TRUE.
2111	outflow_s = .FALSE.
2112	ELSEIF ( nest_domain ) THEN
2113	nest_bound_s = .TRUE.
2114	nest_bound_n = .TRUE.
2115	ENDIF
2116
2117	DEALLOCATE( f2_sub )
2118
2119	ELSE
2120
2121	CALL restrict_fast( f2, r )
2122
2123	ind_even_odd = even_odd_level(grid_level) ! must be after restrict
2124
2125	ENDIF
2126
2127	p2 = 0.0_wp
2128
2129	!
2130	!-- Repeat the same procedure till the coarsest grid is reached
2131	CALL next_mg_level_fast( f2, p2, p3, r )
2132
2133	ENDIF
2134
2135	!
2136	!-- Now follows the prolongation
2137	IF ( grid_level >= 2 ) THEN
2138
2139	!
2140	!-- Prolongation of the new residual. The values are transferred
2141	!-- from the coarse to the next finer grid.
2142	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
2143
2144	#if defined( __parallel )
2145	!
2146	!-- At this level, the new residual first has to be scattered from
2147	!-- PE0 to the other PEs
2148	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, &
2149	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
2150	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
2151
2152	CALL mg_scatter_fast( p2, p2_sub )
2153
2154	!
2155	!-- Therefore, indices of the previous level have to be changed to
2156	!-- subdomain values in between (otherwise, the prolong routine would
2157	!-- expect the gathered array)
2158	nxl_mg_save = nxl_mg(grid_level-1)
2159	nxr_mg_save = nxr_mg(grid_level-1)
2160	nys_mg_save = nys_mg(grid_level-1)
2161	nyn_mg_save = nyn_mg(grid_level-1)
2162	nzt_mg_save = nzt_mg(grid_level-1)
2163	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
2164	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
2165	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
2166	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
2167	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
2168
2169	!
2170	!-- Set switch for routine exchange_horiz, that ghostpoint exchange
2171	!-- has to be carried again out from now on
2172	mg_switch_to_pe0 = .FALSE.
2173
2174	!
2175	!-- For non-cyclic lateral boundary conditions, restore the
2176	!-- in-/outflow conditions
2177	inflow_l = .FALSE.; inflow_r = .FALSE.
2178	inflow_n = .FALSE.; inflow_s = .FALSE.
2179	outflow_l = .FALSE.; outflow_r = .FALSE.
2180	outflow_n = .FALSE.; outflow_s = .FALSE.
2181
2182	IF ( pleft == MPI_PROC_NULL ) THEN
2183	IF ( bc_lr_dirrad ) THEN
2184	inflow_l = .TRUE.
2185	ELSEIF ( bc_lr_raddir ) THEN
2186	outflow_l = .TRUE.
2187	ELSEIF ( nest_domain ) THEN
2188	nest_bound_l = .TRUE.
2189	ENDIF
2190	ENDIF
2191
2192	IF ( pright == MPI_PROC_NULL ) THEN
2193	IF ( bc_lr_dirrad ) THEN
2194	outflow_r = .TRUE.
2195	ELSEIF ( bc_lr_raddir ) THEN
2196	inflow_r = .TRUE.
2197	ELSEIF ( nest_domain ) THEN
2198	nest_bound_r = .TRUE.
2199	ENDIF
2200	ENDIF
2201
2202	IF ( psouth == MPI_PROC_NULL ) THEN
2203	IF ( bc_ns_dirrad ) THEN
2204	outflow_s = .TRUE.
2205	ELSEIF ( bc_ns_raddir ) THEN
2206	inflow_s = .TRUE.
2207	ELSEIF ( nest_domain ) THEN
2208	nest_bound_s = .TRUE.
2209	ENDIF
2210	ENDIF
2211
2212	IF ( pnorth == MPI_PROC_NULL ) THEN
2213	IF ( bc_ns_dirrad ) THEN
2214	inflow_n = .TRUE.
2215	ELSEIF ( bc_ns_raddir ) THEN
2216	outflow_n = .TRUE.
2217	ELSEIF ( nest_domain ) THEN
2218	nest_bound_n = .TRUE.
2219	ENDIF
2220	ENDIF
2221
2222	CALL prolong_fast( p2_sub, p3 )
2223
2224	!
2225	!-- Restore the correct indices of the previous level
2226	nxl_mg(grid_level-1) = nxl_mg_save
2227	nxr_mg(grid_level-1) = nxr_mg_save
2228	nys_mg(grid_level-1) = nys_mg_save
2229	nyn_mg(grid_level-1) = nyn_mg_save
2230	nzt_mg(grid_level-1) = nzt_mg_save
2231
2232	DEALLOCATE( p2_sub )
2233	#endif
2234
2235	ELSE
2236
2237	CALL prolong_fast( p2, p3 )
2238
2239	ENDIF
2240
2241	!
2242	!-- Computation of the new pressure correction. Therefore,
2243	!-- values from prior grids are added up automatically stage by stage.
2244	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
2245	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
2246	DO k = nzb, nzt_mg(grid_level)+1
2247	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
2248	ENDDO
2249	ENDDO
2250	ENDDO
2251
2252	!
2253	!-- Relaxation of the new solution
2254	CALL redblack_fast( f_mg, p_mg )
2255
2256	ENDIF
2257
2258
2259	!
2260	!-- The following few lines serve the steering of the multigrid scheme
2261	IF ( grid_level == maximum_grid_level ) THEN
2262
2263	GOTO 20
2264
2265	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
2266	grid_level_count(grid_level) /= gamma_mg ) THEN
2267
2268	GOTO 10
2269
2270	ENDIF
2271
2272	!
2273	!-- Reset counter for the next call of poismg_fast
2274	grid_level_count(grid_level) = 0
2275
2276	!
2277	!-- Continue with the next finer level. nxl..nzt have to be
2278	!-- set to the finer grid values, because these variables are needed for the
2279	!-- exchange of ghost points in routine exchange_horiz
2280	grid_level = grid_level + 1
2281	ind_even_odd = even_odd_level(grid_level)
2282
2283	nxl = nxl_mg(grid_level)
2284	nxr = nxr_mg(grid_level)
2285	nys = nys_mg(grid_level)
2286	nyn = nyn_mg(grid_level)
2287	nzt = nzt_mg(grid_level)
2288
2289	20 CONTINUE
2290
2291	END SUBROUTINE next_mg_level_fast
2292
2293
2294	!------------------------------------------------------------------------------!
2295	! Description:
2296	! ------------
2297	!> Initial settings for sorting k-dimension from sequential order (alternate
2298	!> even/odd) into blocks of even and odd or vice versa
2299	!------------------------------------------------------------------------------!
2300	SUBROUTINE init_even_odd_blocks
2301
2302
2303	USE arrays_3d, &
2304	ONLY: f1_mg, f2_mg, f3_mg
2305
2306	USE control_parameters, &
2307	ONLY: grid_level, masking_method, maximum_grid_level, topography
2308
2309	USE indices, &
2310	ONLY: nzb, nzt, nzt_mg
2311
2312	USE indices, &
2313	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, &
2314	wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7, &
2315	wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, &
2316	nys_mg, nyn_mg, nzb, nzt_mg
2317
2318	IMPLICIT NONE
2319	!
2320	!-- Local variables
2321	INTEGER(iwp) :: i !<
2322	INTEGER(iwp) :: l !<
2323
2324	LOGICAL, SAVE :: lfirst = .TRUE.
2325
2326
2327	IF ( .NOT. lfirst ) RETURN
2328
2329	ALLOCATE( even_odd_level(maximum_grid_level) )
2330
2331	ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level), &
2332	f2_mg_b(nzb:nzt+1,maximum_grid_level), &
2333	f3_mg_b(nzb:nzt+1,maximum_grid_level) )
2334
2335	!
2336	!-- Set border index between the even and odd block
2337	DO i = maximum_grid_level, 1, -1
2338	even_odd_level(i) = nzt_mg(i) / 2
2339	ENDDO
2340
2341	!
2342	!-- Sort grid coefficients used in red/black scheme and for calculating the
2343	!-- residual to block (even/odd) structure
2344	DO l = maximum_grid_level, 1 , -1
2345	CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l), &
2346	f1_mg_b(nzb:nzt_mg(grid_level)+1,l), &
2347	l )
2348	CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l), &
2349	f2_mg_b(nzb:nzt_mg(grid_level)+1,l), &
2350	l )
2351	CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l), &
2352	f3_mg_b(nzb:nzt_mg(grid_level)+1,l), &
2353	l )
2354	ENDDO
2355
2356	!
2357	!-- Sort flags array for all levels to block (even/odd) structure
2358	IF ( topography /= 'flat' .AND. .NOT. masking_method ) THEN
2359
2360	DO l = maximum_grid_level, 1 , -1
2361	!
2362	!-- Assign the flag level to be calculated
2363	SELECT CASE ( l )
2364	CASE ( 1 )
2365	flags => wall_flags_1
2366	CASE ( 2 )
2367	flags => wall_flags_2
2368	CASE ( 3 )
2369	flags => wall_flags_3
2370	CASE ( 4 )
2371	flags => wall_flags_4
2372	CASE ( 5 )
2373	flags => wall_flags_5
2374	CASE ( 6 )
2375	flags => wall_flags_6
2376	CASE ( 7 )
2377	flags => wall_flags_7
2378	CASE ( 8 )
2379	flags => wall_flags_8
2380	CASE ( 9 )
2381	flags => wall_flags_9
2382	CASE ( 10 )
2383	flags => wall_flags_10
2384	END SELECT
2385
2386	CALL sort_k_to_even_odd_blocks( flags, l )
2387
2388	ENDDO
2389
2390	ENDIF
2391
2392	lfirst = .FALSE.
2393
2394	END SUBROUTINE init_even_odd_blocks
2395
2396
2397	!------------------------------------------------------------------------------!
2398	! Description:
2399	! ------------
2400	!> Special exchange_horiz subroutine for use in redblack. Transfers only
2401	!> "red" or "black" data points.
2402	!------------------------------------------------------------------------------!
2403	SUBROUTINE special_exchange_horiz ( p_mg, color )
2404
2405
2406	USE control_parameters, &
2407	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, &
2408	inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level, &
2409	outflow_l, outflow_n, outflow_r, outflow_s, &
2410	synchronous_exchange
2411
2412	USE indices, &
2413	ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, &
2414	nyn_mg, nzb, nzt, nzt_mg
2415
2416	IMPLICIT NONE
2417
2418	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
2419	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
2420	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
2421
2422	INTEGER(iwp), intent(IN) :: color !<
2423	!
2424	!-- Local variables
2425	INTEGER(iwp) :: i,i1,i2 !<
2426	INTEGER(iwp) :: j,j1,j2 !<
2427	INTEGER(iwp) :: k !<
2428	INTEGER(iwp) :: l !<
2429	INTEGER(iwp) :: jys !<
2430	INTEGER(iwp) :: jyn !<
2431	INTEGER(iwp) :: ixl !<
2432	INTEGER(iwp) :: ixr !<
2433	logical :: sendrecv_in_background_save !<
2434	logical :: synchronous_exchange_save !<
2435
2436	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
2437	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
2438	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: temp !<
2439
2440	#if defined ( __parallel )
2441	sendrecv_in_background_save = sendrecv_in_background
2442	sendrecv_in_background = .FALSE.
2443	synchronous_exchange_save = synchronous_exchange
2444	synchronous_exchange = .FALSE.
2445
2446	l = grid_level
2447
2448	ind_even_odd = even_odd_level(grid_level)
2449
2450	jys = nys_mg(grid_level-1)
2451	jyn = nyn_mg(grid_level-1)
2452	ixl = nxl_mg(grid_level-1)
2453	ixr = nxr_mg(grid_level-1)
2454
2455	!
2456	!-- Restricted transfer only on finer levels with enough data
2457	IF ( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 ) THEN
2458	!
2459	!-- Handling the even k Values
2460	!-- Collecting data for the north - south exchange
2461	!-- Since only every second value has to be transfered, data are stored
2462	!-- on the next coarser grid level, because the arrays on that level
2463	!-- have just the required size
2464	i1 = nxl_mg(grid_level-1)
2465	i2 = nxl_mg(grid_level-1)
2466
2467	DO i = nxl_mg(l), nxr_mg(l), 2
2468	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2469
2470	IF ( j == nys_mg(l) ) THEN
2471	!DIR$ IVDEP
2472	DO k = ind_even_odd+1, nzt_mg(l)
2473	temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
2474	ENDDO
2475	i1 = i1 + 1
2476
2477	ENDIF
2478
2479	IF ( j == nyn_mg(l) ) THEN
2480	!DIR$ IVDEP
2481	DO k = ind_even_odd+1, nzt_mg(l)
2482	temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
2483	ENDDO
2484	i2 = i2 + 1
2485
2486	ENDIF
2487
2488	ENDDO
2489	ENDDO
2490
2491	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2492	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2493
2494	IF ( j == nys_mg(l) ) THEN
2495	!DIR$ IVDEP
2496	DO k = ind_even_odd+1, nzt_mg(l)
2497	temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
2498	ENDDO
2499	i1 = i1 + 1
2500
2501	ENDIF
2502
2503	IF ( j == nyn_mg(l) ) THEN
2504	!DIR$ IVDEP
2505	DO k = ind_even_odd+1, nzt_mg(l)
2506	temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
2507	ENDDO
2508	i2 = i2 + 1
2509
2510	ENDIF
2511
2512	ENDDO
2513	ENDDO
2514
2515	grid_level = grid_level-1
2516
2517	nxl = nxl_mg(grid_level)
2518	nys = nys_mg(grid_level)
2519	nxr = nxr_mg(grid_level)
2520	nyn = nyn_mg(grid_level)
2521	nzt = nzt_mg(grid_level)
2522
2523	send_receive = 'ns'
2524	CALL exchange_horiz( temp, 1 )
2525
2526	grid_level = grid_level+1
2527
2528	i1 = nxl_mg(grid_level-1)
2529	i2 = nxl_mg(grid_level-1)
2530
2531	DO i = nxl_mg(l), nxr_mg(l), 2
2532	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2533
2534	IF ( j == nys_mg(l) ) THEN
2535	!DIR$ IVDEP
2536	DO k = ind_even_odd+1, nzt_mg(l)
2537	p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
2538	ENDDO
2539	i1 = i1 + 1
2540
2541	ENDIF
2542
2543	IF ( j == nyn_mg(l) ) THEN
2544	!DIR$ IVDEP
2545	DO k = ind_even_odd+1, nzt_mg(l)
2546	p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
2547	ENDDO
2548	i2 = i2 + 1
2549
2550	ENDIF
2551
2552	ENDDO
2553	ENDDO
2554
2555	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2556	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2557
2558	IF ( j == nys_mg(l) ) THEN
2559	!DIR$ IVDEP
2560	DO k = ind_even_odd+1, nzt_mg(l)
2561	p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
2562	ENDDO
2563	i1 = i1 + 1
2564
2565	ENDIF
2566
2567	IF ( j == nyn_mg(l) ) THEN
2568	!DIR$ IVDEP
2569	DO k = ind_even_odd+1, nzt_mg(l)
2570	p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
2571	ENDDO
2572	i2 = i2 + 1
2573
2574	ENDIF
2575
2576	ENDDO
2577	ENDDO
2578
2579	!
2580	!-- Collecting data for the left - right exchange
2581	!-- Since only every second value has to be transfered, data are stored
2582	!-- on the next coarser grid level, because the arrays on that level
2583	!-- have just the required size
2584	j1 = nys_mg(grid_level-1)
2585	j2 = nys_mg(grid_level-1)
2586
2587	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2588	DO i = nxl_mg(l), nxr_mg(l), 2
2589
2590	IF ( i == nxl_mg(l) ) THEN
2591	!DIR$ IVDEP
2592	DO k = ind_even_odd+1, nzt_mg(l)
2593	temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
2594	ENDDO
2595	j1 = j1 + 1
2596
2597	ENDIF
2598
2599	IF ( i == nxr_mg(l) ) THEN
2600	!DIR$ IVDEP
2601	DO k = ind_even_odd+1, nzt_mg(l)
2602	temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
2603	ENDDO
2604	j2 = j2 + 1
2605
2606	ENDIF
2607
2608	ENDDO
2609	ENDDO
2610
2611	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2612	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2613
2614	IF ( i == nxl_mg(l) ) THEN
2615	!DIR$ IVDEP
2616	DO k = ind_even_odd+1, nzt_mg(l)
2617	temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
2618	ENDDO
2619	j1 = j1 + 1
2620
2621	ENDIF
2622
2623	IF ( i == nxr_mg(l) ) THEN
2624	!DIR$ IVDEP
2625	DO k = ind_even_odd+1, nzt_mg(l)
2626	temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
2627	ENDDO
2628	j2 = j2 + 1
2629
2630	ENDIF
2631
2632	ENDDO
2633	ENDDO
2634
2635	grid_level = grid_level-1
2636	send_receive = 'lr'
2637
2638	CALL exchange_horiz( temp, 1 )
2639
2640	grid_level = grid_level+1
2641
2642	j1 = nys_mg(grid_level-1)
2643	j2 = nys_mg(grid_level-1)
2644
2645	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2646	DO i = nxl_mg(l), nxr_mg(l), 2
2647
2648	IF ( i == nxl_mg(l) ) THEN
2649	!DIR$ IVDEP
2650	DO k = ind_even_odd+1, nzt_mg(l)
2651	p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1)
2652	ENDDO
2653	j1 = j1 + 1
2654
2655	ENDIF
2656
2657	IF ( i == nxr_mg(l) ) THEN
2658	!DIR$ IVDEP
2659	DO k = ind_even_odd+1, nzt_mg(l)
2660	p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1)
2661	ENDDO
2662	j2 = j2 + 1
2663
2664	ENDIF
2665
2666	ENDDO
2667	ENDDO
2668
2669	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2670	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2671
2672	IF ( i == nxl_mg(l) ) THEN
2673	!DIR$ IVDEP
2674	DO k = ind_even_odd+1, nzt_mg(l)
2675	p_mg(k,j,nxr_mg(l)+1) = temp(k-ind_even_odd,j1,ixr+1)
2676	ENDDO
2677	j1 = j1 + 1
2678
2679	ENDIF
2680
2681	IF ( i == nxr_mg(l) ) THEN
2682	!DIR$ IVDEP
2683	DO k = ind_even_odd+1, nzt_mg(l)
2684	p_mg(k,j,nxl_mg(l)-1) = temp(k-ind_even_odd,j2,ixl-1)
2685	ENDDO
2686	j2 = j2 + 1
2687
2688	ENDIF
2689
2690	ENDDO
2691	ENDDO
2692
2693	!
2694	!-- Now handling the even k values
2695	!-- Collecting data for the north - south exchange
2696	!-- Since only every second value has to be transfered, data are stored
2697	!-- on the next coarser grid level, because the arrays on that level
2698	!-- have just the required size
2699	i1 = nxl_mg(grid_level-1)
2700	i2 = nxl_mg(grid_level-1)
2701
2702	DO i = nxl_mg(l), nxr_mg(l), 2
2703	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2704
2705	IF ( j == nys_mg(l) ) THEN
2706	!DIR$ IVDEP
2707	DO k = nzb+1, ind_even_odd
2708	temp(k,jys,i1) = p_mg(k,j,i)
2709	ENDDO
2710	i1 = i1 + 1
2711
2712	ENDIF
2713
2714	IF ( j == nyn_mg(l) ) THEN
2715	!DIR$ IVDEP
2716	DO k = nzb+1, ind_even_odd
2717	temp(k,jyn,i2) = p_mg(k,j,i)
2718	ENDDO
2719	i2 = i2 + 1
2720
2721	ENDIF
2722
2723	ENDDO
2724	ENDDO
2725
2726	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2727	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2728
2729	IF ( j == nys_mg(l) ) THEN
2730	!DIR$ IVDEP
2731	DO k = nzb+1, ind_even_odd
2732	temp(k,jys,i1) = p_mg(k,j,i)
2733	ENDDO
2734	i1 = i1 + 1
2735
2736	ENDIF
2737
2738	IF ( j == nyn_mg(l) ) THEN
2739	!DIR$ IVDEP
2740	DO k = nzb+1, ind_even_odd
2741	temp(k,jyn,i2) = p_mg(k,j,i)
2742	ENDDO
2743	i2 = i2 + 1
2744
2745	ENDIF
2746
2747	ENDDO
2748	ENDDO
2749
2750	grid_level = grid_level-1
2751
2752	send_receive = 'ns'
2753	CALL exchange_horiz( temp, 1 )
2754
2755	grid_level = grid_level+1
2756
2757	i1 = nxl_mg(grid_level-1)
2758	i2 = nxl_mg(grid_level-1)
2759
2760	DO i = nxl_mg(l), nxr_mg(l), 2
2761	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2762
2763	IF ( j == nys_mg(l) ) THEN
2764	!DIR$ IVDEP
2765	DO k = nzb+1, ind_even_odd
2766	p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
2767	ENDDO
2768	i1 = i1 + 1
2769
2770	ENDIF
2771
2772	IF ( j == nyn_mg(l) ) THEN
2773	!DIR$ IVDEP
2774	DO k = nzb+1, ind_even_odd
2775	p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
2776	ENDDO
2777	i2 = i2 + 1
2778
2779	ENDIF
2780
2781	ENDDO
2782	ENDDO
2783
2784	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2785	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2786
2787	IF ( j == nys_mg(l) ) THEN
2788	!DIR$ IVDEP
2789	DO k = nzb+1, ind_even_odd
2790	p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
2791	ENDDO
2792	i1 = i1 + 1
2793
2794	ENDIF
2795
2796	IF ( j == nyn_mg(l) ) THEN
2797	!DIR$ IVDEP
2798	DO k = nzb+1, ind_even_odd
2799	p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
2800	ENDDO
2801	i2 = i2 + 1
2802
2803	ENDIF
2804
2805	ENDDO
2806	ENDDO
2807
2808	j1 = nys_mg(grid_level-1)
2809	j2 = nys_mg(grid_level-1)
2810
2811	DO i = nxl_mg(l), nxr_mg(l), 2
2812	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2813
2814	IF ( i == nxl_mg(l) ) THEN
2815	!DIR$ IVDEP
2816	DO k = nzb+1, ind_even_odd
2817	temp(k,j1,ixl) = p_mg(k,j,i)
2818	ENDDO
2819	j1 = j1 + 1
2820
2821	ENDIF
2822
2823	IF ( i == nxr_mg(l) ) THEN
2824	!DIR$ IVDEP
2825	DO k = nzb+1, ind_even_odd
2826	temp(k,j2,ixr) = p_mg(k,j,i)
2827	ENDDO
2828	j2 = j2 + 1
2829
2830	ENDIF
2831
2832	ENDDO
2833	ENDDO
2834
2835	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2836	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2837
2838	IF ( i == nxl_mg(l) ) THEN
2839	!DIR$ IVDEP
2840	DO k = nzb+1, ind_even_odd
2841	temp(k,j1,ixl) = p_mg(k,j,i)
2842	ENDDO
2843	j1 = j1 + 1
2844
2845	ENDIF
2846
2847	IF ( i == nxr_mg(l) ) THEN
2848	!DIR$ IVDEP
2849	DO k = nzb+1, ind_even_odd
2850	temp(k,j2,ixr) = p_mg(k,j,i)
2851	ENDDO
2852	j2 = j2 + 1
2853
2854	ENDIF
2855
2856	ENDDO
2857	ENDDO
2858
2859	grid_level = grid_level-1
2860
2861	send_receive = 'lr'
2862	CALL exchange_horiz( temp, 1 )
2863
2864	grid_level = grid_level+1
2865
2866	nxl = nxl_mg(grid_level)
2867	nys = nys_mg(grid_level)
2868	nxr = nxr_mg(grid_level)
2869	nyn = nyn_mg(grid_level)
2870	nzt = nzt_mg(grid_level)
2871
2872	j1 = nys_mg(grid_level-1)
2873	j2 = nys_mg(grid_level-1)
2874
2875	DO i = nxl_mg(l), nxr_mg(l), 2
2876	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2877
2878	IF ( i == nxl_mg(l) ) THEN
2879	!DIR$ IVDEP
2880	DO k = nzb+1, ind_even_odd
2881	p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1)
2882	ENDDO
2883	j1 = j1 + 1
2884
2885	ENDIF
2886
2887	IF ( i == nxr_mg(l) ) THEN
2888	!DIR$ IVDEP
2889	DO k = nzb+1, ind_even_odd
2890	p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1)
2891	ENDDO
2892	j2 = j2 + 1
2893
2894	ENDIF
2895
2896	ENDDO
2897	ENDDO
2898
2899	DO i = nxl_mg(l)+1, nxr_mg(l), 2
2900	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2901
2902	IF ( i == nxl_mg(l) ) THEN
2903	!DIR$ IVDEP
2904	DO k = nzb+1, ind_even_odd
2905	p_mg(k,j,nxr_mg(l)+1) = temp(k,j1,ixr+1)
2906	ENDDO
2907	j1 = j1 + 1
2908
2909	ENDIF
2910
2911	IF ( i == nxr_mg(l) ) THEN
2912	!DIR$ IVDEP
2913	DO k = nzb+1, ind_even_odd
2914	p_mg(k,j,nxl_mg(l)-1) = temp(k,j2,ixl-1)
2915	ENDDO
2916	j2 = j2 + 1
2917
2918	ENDIF
2919
2920	ENDDO
2921	ENDDO
2922
2923	ELSE
2924
2925	!
2926	!-- Standard horizontal ghost boundary exchange for small coarse grid
2927	!-- levels, where the transfer time is latency bound
2928	CALL exchange_horiz( p_mg, 1 )
2929
2930	ENDIF
2931
2932	!
2933	!-- Reset values to default PALM setup
2934	sendrecv_in_background = sendrecv_in_background_save
2935	synchronous_exchange = synchronous_exchange_save
2936	send_receive = 'al'
2937	#else
2938
2939	!
2940	!-- Standard horizontal ghost boundary exchange for small coarse grid
2941	!-- levels, where the transfer time is latency bound
2942	CALL exchange_horiz( p_mg, 1 )
2943	#endif
2944
2945	END SUBROUTINE special_exchange_horiz
2946
2947	END MODULE poismg_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |