Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |