Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |