Home

Context Navigation

source: palm/trunk/SOURCE/poismg_fast_mod.f90 @ 1883

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |