Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |