Home

Context Navigation

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

Last change on this file since 3574 was 3534, checked in by raasch, 6 years ago
inifor integrated in build mechanism, some bugfixes in inifor to avoid compile time errors, batch_scp for sending back the job protocol file is called via login-node if a login-node has been set in the config-file, ssh-calls rearranged to avoid output of system/user-profile messages
Property svn:keywords set to `Id`
File size: 86.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |