Home

Context Navigation

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

Last change on this file since 1754 was 1683, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 114.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |