Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |