Home

Context Navigation

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

Last change on this file since 2941 was 2939, checked in by suehring, 7 years ago
Set lateral boundary conditions for divergence in multigrid solver
Property svn:keywords set to `Id`
File size: 89.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |