Home

Context Navigation

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

Last change on this file since 1883 was 1851, checked in by maronga, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 115.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |