Home

Context Navigation

poismg_mod.f90 @ 4865

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

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

Context Navigation

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

Download in other formats: