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

Rev	Line
[1682]	1	!> @file poismg_fast.f90
[1575]	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	! -----------------
[1610]	21	!
[1683]	22	!
[1576]	23	! Former revisions:
	24	! -----------------
	25	! $Id: poismg_fast.f90 1683 2015-10-07 23:57:51Z maronga $
	26	!
[1683]	27	! 1682 2015-10-07 23:56:08Z knoop
	28	! Code annotations made doxygen readable
	29	!
[1610]	30	! 1609 2015-07-03 15:37:58Z maronga
	31	! Bugfix: allow compilation without __parallel.
	32	!
[1576]	33	! 1575 2015-03-27 09:56:27Z raasch
[1575]	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	! ------------
[1682]	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.
[1575]	60	!------------------------------------------------------------------------------!
[1682]	61	MODULE poismg_mod
	62
[1575]	63
	64	USE cpulog, &
	65	ONLY: cpu_log, log_point_s
[1609]	66
[1575]	67	USE kinds
[1609]	68
[1575]	69	USE pegrid
	70
	71	PRIVATE
	72
[1682]	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
[1575]	75
[1682]	76	REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE :: f1_mg_b, f2_mg_b, f3_mg_b !< blocked version of f1_mg ...
[1575]	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
[1682]	92	!------------------------------------------------------------------------------!
	93	! Description:
	94	! ------------
	95	!> Solves the Poisson equation for the perturbation pressure with a multigrid
	96	!> V- or W-Cycle scheme.
	97	!------------------------------------------------------------------------------!
[1575]	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
[1682]	117	REAL(wp) :: maxerror !<
	118	REAL(wp) :: maximum_mgcycles !<
	119	REAL(wp) :: residual_norm !<
[1575]	120
[1682]	121	REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !<
[1575]	122
[1682]	123	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !<
[1575]	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	! ------------
[1682]	232	!> Computes the residual of the perturbation pressure.
[1575]	233	!------------------------------------------------------------------------------!
[1682]	234	SUBROUTINE resid_fast( f_mg, p_mg, r )
[1575]	235
[1682]	236
[1575]	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
[1682]	260	INTEGER(iwp) :: km1 !<
	261	INTEGER(iwp) :: kp1 !<
[1575]	262
	263	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	264	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	265	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1575]	266	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	267	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	268	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1575]	269	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	270	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	271	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
[1575]	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	! ------------
[1682]	428	!> Interpolates the residual on the next coarser grid with "full weighting"
	429	!> scheme
[1575]	430	!------------------------------------------------------------------------------!
[1682]	431	SUBROUTINE restrict_fast( f_mg, r )
[1575]	432
[1682]	433
[1575]	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
[1682]	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 !<
[1575]	456
[1682]	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 !<
[1575]	474
	475	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	476	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	477	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1575]	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, &
[1682]	481	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !<
[1575]	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	! ------------
[1682]	673	!> Interpolates the correction of the perturbation pressure
	674	!> to the next finer grid.
[1575]	675	!------------------------------------------------------------------------------!
[1682]	676	SUBROUTINE prolong_fast( p, temp )
[1575]	677
[1682]	678
[1575]	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
[1682]	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
[1575]	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, &
[1682]	699	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !<
[1575]	700
	701	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	702	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	703	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !<
[1575]	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	! ------------
[1682]	823	!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
	824	!> 3D-Red-Black decomposition (GS-RB) is used.
[1575]	825	!------------------------------------------------------------------------------!
[1682]	826	SUBROUTINE redblack_fast( f_mg, p_mg )
[1575]	827
[1682]	828
[1575]	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
[1682]	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 !<
[1575]	859
[1682]	860	LOGICAL :: unroll !<
[1575]	861
[1682]	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 !<
[1575]	868
	869	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	870	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	871	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1575]	872	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	873	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	874	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1575]	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	! ------------
[1682]	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
[1575]	1536	!------------------------------------------------------------------------------!
[1682]	1537	SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel )
[1575]	1538
[1682]	1539
[1575]	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, &
[1682]	1552	nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: p_mg !<
[1575]	1553	!
	1554	!-- Local variables
[1682]	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 !<
[1575]	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	! ------------
[1682]	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
[1575]	1604	!------------------------------------------------------------------------------!
[1682]	1605	SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel )
[1575]	1606
[1682]	1607
[1575]	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
[1682]	1620	INTEGER(iwp) :: ind !<
	1621	INTEGER(iwp) :: k !<
[1575]	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	! ------------
[1682]	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
[1575]	1651	!------------------------------------------------------------------------------!
[1682]	1652	SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
[1575]	1653
[1682]	1654
[1575]	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, &
[1682]	1667	nxl_mg(glevel)-1:nxr_mg(glevel)+1) :: i_mg !<
[1575]	1668	!
	1669	!-- Local variables
[1682]	1670	INTEGER(iwp) :: i !<
	1671	INTEGER(iwp) :: j !<
	1672	INTEGER(iwp) :: k !<
	1673	INTEGER(iwp) :: l !<
	1674	INTEGER(iwp) :: ind !<
[1575]	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
[1609]	1696	#if defined ( __parallel )
[1575]	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
[1609]	1701	#endif
[1575]	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	! ------------
[1682]	1722	!> Sort k-dimension from blocks of even and odd into sequential
[1575]	1723	!------------------------------------------------------------------------------!
[1682]	1724	SUBROUTINE sort_k_to_sequential( p_mg )
[1575]	1725
[1682]	1726
[1575]	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, &
[1682]	1737	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1575]	1738	!
	1739	!-- Local variables
[1682]	1740	INTEGER(iwp) :: i !<
	1741	INTEGER(iwp) :: j !<
	1742	INTEGER(iwp) :: k !<
	1743	INTEGER(iwp) :: l !<
	1744	INTEGER(iwp) :: ind !<
[1575]	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
[1682]	1774	!------------------------------------------------------------------------------!
	1775	! Description:
	1776	! ------------
	1777	!> Gather subdomain data from all PEs.
	1778	!------------------------------------------------------------------------------!
[1575]	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
[1682]	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 !<
[1575]	1800
	1801	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1802	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1803	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !<
[1575]	1804	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1805	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1806	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !<
[1575]	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, &
[1682]	1810	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !<
[1575]	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
[1682]	1854	!------------------------------------------------------------------------------!
	1855	! Description:
	1856	! ------------
	1857	!> @todo It might be possible to improve the speed of this routine by using
	1858	!> non-blocking communication
	1859	!------------------------------------------------------------------------------!
[1575]	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
[1682]	1873	INTEGER(iwp) :: nwords !<
[1575]	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, &
[1682]	1877	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[1575]	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, &
[1682]	1881	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !<
[1575]	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	! ------------
[1682]	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.
[1575]	1908	!------------------------------------------------------------------------------!
[1682]	1909	RECURSIVE SUBROUTINE next_mg_level_fast( f_mg, p_mg, p3, r )
[1575]	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
[1682]	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 !<
[1575]	1932
	1933	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1934	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1935	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1575]	1936	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1937	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1938	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1575]	1939	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1940	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1941	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !<
[1575]	1942	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1943	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1944	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
[1575]	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, &
[1682]	1948	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !<
[1575]	1949	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1950	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1951	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[1575]	1952
[1682]	1953	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !<
	1954	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !<
[1575]	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	! ------------
[1682]	2250	!> Initial settings for sorting k-dimension from sequential order (alternate
	2251	!> even/odd) into blocks of even and odd or vice versa
[1575]	2252	!------------------------------------------------------------------------------!
[1682]	2253	SUBROUTINE init_even_odd_blocks
[1575]	2254
[1682]	2255
[1575]	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
[1682]	2274	INTEGER(iwp) :: i !<
	2275	INTEGER(iwp) :: l !<
[1575]	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	! ------------
[1682]	2353	!> Special exchange_horiz subroutine for use in redblack. Transfers only
	2354	!> "red" or "black" data points.
[1575]	2355	!------------------------------------------------------------------------------!
[1682]	2356	SUBROUTINE special_exchange_horiz ( p_mg, color )
[1575]	2357
[1682]	2358
[1575]	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, &
[1682]	2373	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1575]	2374
[1682]	2375	INTEGER(iwp), intent(IN) :: color !<
[1575]	2376	!
	2377	!-- Local variables
[1682]	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 !<
[1575]	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, &
[1682]	2391	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: temp !<
[1575]	2392
[1609]	2393	#if defined ( __parallel )
[1575]	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
[1609]	2877
[1575]	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'
[1609]	2890	#else
[1575]	2891
[1609]	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
[1575]	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 |