Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |