Home

Context Navigation

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

Last change on this file since 2005 was 2001, checked in by knoop, 9 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 86.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |