Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |