Home

Context Navigation

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

Last change on this file since 4680 was 4649, checked in by raasch, 4 years ago
files re-formatted to follow the PALM coding standard
Property svn:keywords set to `Id`
File size: 94.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |