Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |