Home

Context Navigation

source: palm/trunk/SOURCE/poismg_mod.f90 @ 1941

Last change on this file since 1941 was 1935, checked in by hellstea, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |