Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |