Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |