Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |