Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |