Home

Context Navigation

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

Last change on this file since 2849 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 88.4 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |