Home

Context Navigation

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

Last change on this file since 2312 was 2298, checked in by raasch, 7 years ago
write_binary is of type LOGICAL now, MPI2-related code removed, obsolete variables removed, sendrecv_in_background related parts removed, missing variable descriptions added
Property svn:keywords set to `Id`
File size: 87.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |