Home

Context Navigation

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

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |