Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |