Home

Context Navigation

source: palm/trunk/SOURCE/poismg_fast.f90 @ 1788

Last change on this file since 1788 was 1763, checked in by hellstea, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 115.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |