Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |