Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |