Home

Context Navigation

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

Last change on this file since 2186 was 2101, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 88.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |