Home

Context Navigation

source: palm/trunk/SOURCE/poismg.f90 @ 1682

Last change on this file since 1682 was 1682, checked in by knoop, 9 years ago
Code annotations made doxygen readable
Property svn:keywords set to `Id`
File size: 65.2 KB

Rev	Line
[1682]	1	!> @file poismg.f90
[1036]	2	!--------------------------------------------------------------------------------!
	3	! This file is part of PALM.
	4	!
	5	! PALM is free software: you can redistribute it and/or modify it under the terms
	6	! of the GNU General Public License as published by the Free Software Foundation,
	7	! either version 3 of the License, or (at your option) any later version.
	8	!
	9	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
	10	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	11	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
	12	!
	13	! You should have received a copy of the GNU General Public License along with
	14	! PALM. If not, see <http://www.gnu.org/licenses/>.
	15	!
[1310]	16	! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]	17	!--------------------------------------------------------------------------------!
	18	!
[257]	19	! Current revisions:
[1]	20	! -----------------
[1682]	21	! Code annotations made doxygen readable
[1354]	22	!
[1321]	23	! Former revisions:
	24	! -----------------
	25	! $Id: poismg.f90 1682 2015-10-07 23:56:08Z knoop $
	26	!
[1354]	27	! 1353 2014-04-08 15:21:23Z heinze
	28	! REAL constants provided with KIND-attribute
	29	!
[1323]	30	! 1322 2014-03-20 16:38:49Z raasch
	31	! REAL constants defined as wp-kind
	32	!
[1321]	33	! 1320 2014-03-20 08:40:49Z raasch
[1320]	34	! ONLY-attribute added to USE-statements,
	35	! kind-parameters added to all INTEGER and REAL declaration statements,
	36	! kinds are defined in new module kinds,
	37	! old module precision_kind is removed,
	38	! revision history before 2012 removed,
	39	! comment fields (!:) to be used for variable explanations added to
	40	! all variable declaration statements
[708]	41	!
[1319]	42	! 1318 2014-03-17 13:35:16Z raasch
	43	! module interfaces removed
	44	!
[1160]	45	! 1159 2013-05-21 11:58:22Z fricke
	46	! bc_lr/ns_dirneu/neudir removed
	47	!
[1093]	48	! 1092 2013-02-02 11:24:22Z raasch
	49	! unused variables removed
	50	!
[1057]	51	! 1056 2012-11-16 15:28:04Z raasch
	52	! Bugfix: all ghost points have to be used for allocating p3
	53	! arrays p2, f2, and f2_l changed from allocatable to automatic
	54	!
[1037]	55	! 1036 2012-10-22 13:43:42Z raasch
	56	! code put under GPL (PALM 3.9)
	57	!
[997]	58	! 996 2012-09-07 10:41:47Z raasch
	59	! little reformatting
	60	!
[979]	61	! 978 2012-08-09 08:28:32Z fricke
	62	! bc_lr/ns_dirneu/neudir added
	63	!
[881]	64	! 880 2012-04-13 06:28:59Z raasch
	65	! Bugfix: preprocessor statements for parallel execution added
	66	!
[779]	67	! 778 2011-11-07 14:18:25Z fricke
	68	! Allocation of p3 changes when multigrid is used and the collected field on PE0
	69	! has more grid points than the subdomain of an PE.
	70	!
[708]	71	! 707 2011-03-29 11:39:40Z raasch
[707]	72	! p_loc is used instead of p in the main routine (poismg).
	73	! On coarse grid levels, gathered data are identically processed on all PEs
	74	! (before, on PE0 only), so that the subsequent scattering of data is not
	75	! neccessary any more.
	76	! bc_lr/ns replaced by bc_lr/ns_cyc/dirrad/raddir
	77	! Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
	78	! resid and restrict. They were missed before which may have led to
	79	! unpredictable results.
[1]	80	!
[668]	81	! 667 2010-12-23 12:06:00Z suehring/gryschka
	82	! Calls of exchange_horiz are modified.
	83	!
[623]	84	! 622 2010-12-10 08:08:13Z raasch
	85	! optional barriers included in order to speed up collective operations
	86	!
[392]	87	! 257 2009-03-11 15:17:42Z heinze
	88	! Output of messages replaced by message handling routine.
	89	!
[198]	90	! 181 2008-07-30 07:07:47Z raasch
	91	! Bugfix: grid_level+1 has to be used in restrict for flags-array
	92	!
[139]	93	! 114 2007-10-10 00:03:15Z raasch
	94	! Boundary conditions at walls are implicitly set using flag arrays. Only
	95	! Neumann BC is allowed. Upper walls are still not realized.
	96	! Bottom and top BCs for array f_mg in restrict removed because boundary
	97	! values are not needed (right hand side of SOR iteration).
	98	!
[77]	99	! 75 2007-03-22 09:54:05Z raasch
	100	! 2nd+3rd argument removed from exchange horiz
	101	!
[3]	102	! RCS Log replace by Id keyword, revision history cleaned up
	103	!
[1]	104	! Revision 1.6 2005/03/26 20:55:54 raasch
	105	! Implementation of non-cyclic (Neumann) horizontal boundary conditions,
	106	! routine prolong simplified (one call of exchange_horiz spared)
	107	!
	108	! Revision 1.1 2001/07/20 13:10:51 raasch
	109	! Initial revision
	110	!
	111	!
	112	! Description:
	113	! ------------
[1682]	114	!> Solves the Poisson equation for the perturbation pressure with a multigrid
	115	!> V- or W-Cycle scheme.
	116	!>
	117	!> This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
	118	!> September 2000 - July 2001.
	119	!>
	120	!> @attention Loop unrolling and cache optimization in SOR-Red/Black method
	121	!> still does not give the expected speedup!
	122	!>
	123	!> @todo Further work required.
[1]	124	!------------------------------------------------------------------------------!
[1682]	125	SUBROUTINE poismg( r )
	126
[1]	127
[1320]	128	USE arrays_3d, &
	129	ONLY: d, p_loc
	130
	131	USE control_parameters, &
	132	ONLY: gathered_size, grid_level, grid_level_count, &
	133	maximum_grid_level, message_string, mgcycles, mg_cycles, &
	134	mg_switch_to_pe0_level, residual_limit, subdomain_size
	135
	136	USE cpulog, &
	137	ONLY: cpu_log, log_point_s
	138
	139	USE indices, &
	140	ONLY: nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn, &
	141	nyng, nyn_mg, nzb, nzt, nzt_mg
	142
	143	USE kinds
	144
[1]	145	USE pegrid
	146
	147	IMPLICIT NONE
	148
[1682]	149	REAL(wp) :: maxerror !<
	150	REAL(wp) :: maximum_mgcycles !<
	151	REAL(wp) :: residual_norm !<
[1]	152
[1682]	153	REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) :: r !<
[1]	154
[1682]	155	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p3 !<
[1]	156
	157
	158	CALL cpu_log( log_point_s(29), 'poismg', 'start' )
	159	!
	160	!-- Initialize arrays and variables used in this subroutine
	161
[879]	162	!-- If the number of grid points of the gathered grid, which is collected
	163	!-- on PE0, is larger than the number of grid points of an PE, than array
	164	!-- p3 will be enlarged.
[778]	165	IF ( gathered_size > subdomain_size ) THEN
[879]	166	ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg( &
	167	mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
	168	nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg( &
[778]	169	mg_switch_to_pe0_level)+1) )
	170	ELSE
[1056]	171	ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[778]	172	ENDIF
[1056]	173
[1353]	174	p3 = 0.0_wp
[879]	175
[1]	176	!
[707]	177	!-- Ghost boundaries have to be added to divergence array.
	178	!-- Exchange routine needs to know the grid level!
	179	grid_level = maximum_grid_level
[667]	180	CALL exchange_horiz( d, 1)
[1]	181	d(nzb,:,:) = d(nzb+1,:,:)
	182
	183	!
	184	!-- Initiation of the multigrid scheme. Does n cycles until the
	185	!-- residual is smaller than the given limit. The accuracy of the solution
	186	!-- of the poisson equation will increase with the number of cycles.
	187	!-- If the number of cycles is preset by the user, this number will be
	188	!-- carried out regardless of the accuracy.
[707]	189	grid_level_count = 0
	190	mgcycles = 0
[1]	191	IF ( mg_cycles == -1 ) THEN
	192	maximum_mgcycles = 0
[1353]	193	residual_norm = 1.0_wp
[1]	194	ELSE
	195	maximum_mgcycles = mg_cycles
[1353]	196	residual_norm = 0.0_wp
[1]	197	ENDIF
	198
	199	DO WHILE ( residual_norm > residual_limit .OR. &
	200	mgcycles < maximum_mgcycles )
[778]	201
	202	CALL next_mg_level( d, p_loc, p3, r)
[1]	203
	204	!
	205	!-- Calculate the residual if the user has not preset the number of
	206	!-- cycles to be performed
	207	IF ( maximum_mgcycles == 0 ) THEN
[707]	208	CALL resid( d, p_loc, r )
[1]	209	maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
[778]	210
[1]	211	#if defined( __parallel )
[622]	212	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[778]	213	CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL, MPI_SUM, &
[1]	214	comm2d, ierr)
	215	#else
[778]	216	residual_norm = maxerror
[1]	217	#endif
	218	residual_norm = SQRT( residual_norm )
	219	ENDIF
	220
	221	mgcycles = mgcycles + 1
	222
	223	!
	224	!-- If the user has not limited the number of cycles, stop the run in case
	225	!-- of insufficient convergence
	226	IF ( mgcycles > 1000 .AND. mg_cycles == -1 ) THEN
[257]	227	message_string = 'no sufficient convergence within 1000 cycles'
	228	CALL message( 'poismg', 'PA0283', 1, 2, 0, 6, 0 )
[1]	229	ENDIF
	230
	231	ENDDO
	232
	233	DEALLOCATE( p3 )
	234
[707]	235	!
	236	!-- Unset the grid level. Variable is used to determine the MPI datatypes for
	237	!-- ghost point exchange
	238	grid_level = 0
	239
[1]	240	CALL cpu_log( log_point_s(29), 'poismg', 'stop' )
	241
	242	END SUBROUTINE poismg
	243
	244
	245	!------------------------------------------------------------------------------!
	246	! Description:
	247	! ------------
[1682]	248	!> Computes the residual of the perturbation pressure.
[1]	249	!------------------------------------------------------------------------------!
[1682]	250	SUBROUTINE resid( f_mg, p_mg, r )
[1]	251
[1682]	252
[1320]	253	USE arrays_3d, &
	254	ONLY: f1_mg, f2_mg, f3_mg
[1]	255
[1320]	256	USE control_parameters, &
	257	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
	258	inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r, &
	259	outflow_s
	260
	261	USE grid_variables, &
	262	ONLY: ddx2_mg, ddy2_mg
	263
	264	USE indices, &
	265	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
	266	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
	267	wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
	268	nzb, nzt_mg
	269
	270	USE kinds
	271
[1]	272	IMPLICIT NONE
	273
[1320]	274	INTEGER(iwp) :: i
	275	INTEGER(iwp) :: j
	276	INTEGER(iwp) :: k
	277	INTEGER(iwp) :: l
[1]	278
[1320]	279	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
[1]	280	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	281	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1320]	282	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	283	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	284	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1320]	285	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	286	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	287	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
[1]	288
	289	!
	290	!-- Calculate the residual
	291	l = grid_level
	292
[114]	293	!
	294	!-- Choose flag array of this level
	295	SELECT CASE ( l )
	296	CASE ( 1 )
	297	flags => wall_flags_1
	298	CASE ( 2 )
	299	flags => wall_flags_2
	300	CASE ( 3 )
	301	flags => wall_flags_3
	302	CASE ( 4 )
	303	flags => wall_flags_4
	304	CASE ( 5 )
	305	flags => wall_flags_5
	306	CASE ( 6 )
	307	flags => wall_flags_6
	308	CASE ( 7 )
	309	flags => wall_flags_7
	310	CASE ( 8 )
	311	flags => wall_flags_8
	312	CASE ( 9 )
	313	flags => wall_flags_9
	314	CASE ( 10 )
	315	flags => wall_flags_10
	316	END SELECT
	317
[1]	318	!$OMP PARALLEL PRIVATE (i,j,k)
	319	!$OMP DO
	320	DO i = nxl_mg(l), nxr_mg(l)
	321	DO j = nys_mg(l), nyn_mg(l)
	322	DO k = nzb+1, nzt_mg(l)
[114]	323	r(k,j,i) = f_mg(k,j,i) &
	324	- ddx2_mg(l) * &
	325	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	326	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	327	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	328	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	329	- ddy2_mg(l) * &
	330	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	331	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	332	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	333	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	334	- f2_mg(k,l) * p_mg(k+1,j,i) &
	335	- f3_mg(k,l) * &
	336	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	337	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
[1]	338	+ f1_mg(k,l) * p_mg(k,j,i)
[114]	339	!
	340	!-- Residual within topography should be zero
[1353]	341	r(k,j,i) = r(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
[1]	342	ENDDO
	343	ENDDO
	344	ENDDO
	345	!$OMP END PARALLEL
	346
	347	!
	348	!-- Horizontal boundary conditions
[667]	349	CALL exchange_horiz( r, 1)
[1]	350
[707]	351	IF ( .NOT. bc_lr_cyc ) THEN
[1]	352	IF ( inflow_l .OR. outflow_l ) r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
	353	IF ( inflow_r .OR. outflow_r ) r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
	354	ENDIF
	355
[707]	356	IF ( .NOT. bc_ns_cyc ) THEN
[1]	357	IF ( inflow_n .OR. outflow_n ) r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
	358	IF ( inflow_s .OR. outflow_s ) r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
	359	ENDIF
	360
	361	!
[707]	362	!-- Boundary conditions at bottom and top of the domain.
	363	!-- These points are not handled by the above loop. Points may be within
	364	!-- buildings, but that doesn't matter.
	365	IF ( ibc_p_b == 1 ) THEN
	366	r(nzb,:,: ) = r(nzb+1,:,:)
	367	ELSE
[1353]	368	r(nzb,:,: ) = 0.0_wp
[707]	369	ENDIF
	370
[1]	371	IF ( ibc_p_t == 1 ) THEN
	372	r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
	373	ELSE
[1353]	374	r(nzt_mg(l)+1,:,: ) = 0.0_wp
[1]	375	ENDIF
	376
	377
	378	END SUBROUTINE resid
	379
	380
	381	!------------------------------------------------------------------------------!
	382	! Description:
	383	! ------------
[1682]	384	!> Interpolates the residual on the next coarser grid with "full weighting"
	385	!> scheme.
[1]	386	!------------------------------------------------------------------------------!
[1682]	387	SUBROUTINE restrict( f_mg, r )
[1]	388
[1682]	389
[1320]	390	USE control_parameters, &
	391	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
	392	inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r, &
	393	outflow_s
[1]	394
[1320]	395	USE indices, &
	396	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
	397	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
	398	wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
	399	nzb, nzt_mg
	400
	401	USE kinds
	402
[1]	403	IMPLICIT NONE
	404
[1682]	405	INTEGER(iwp) :: i !<
	406	INTEGER(iwp) :: ic !<
	407	INTEGER(iwp) :: j !<
	408	INTEGER(iwp) :: jc !<
	409	INTEGER(iwp) :: k !<
	410	INTEGER(iwp) :: kc !<
	411	INTEGER(iwp) :: l !<
[1]	412
[1682]	413	REAL(wp) :: rkjim !<
	414	REAL(wp) :: rkjip !<
	415	REAL(wp) :: rkjmi !<
	416	REAL(wp) :: rkjmim !<
	417	REAL(wp) :: rkjmip !<
	418	REAL(wp) :: rkjpi !<
	419	REAL(wp) :: rkjpim !<
	420	REAL(wp) :: rkjpip !<
	421	REAL(wp) :: rkmji !<
	422	REAL(wp) :: rkmjim !<
	423	REAL(wp) :: rkmjip !<
	424	REAL(wp) :: rkmjmi !<
	425	REAL(wp) :: rkmjmim !<
	426	REAL(wp) :: rkmjmip !<
	427	REAL(wp) :: rkmjpi !<
	428	REAL(wp) :: rkmjpim !<
	429	REAL(wp) :: rkmjpip !<
[114]	430
[1320]	431	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	432	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	433	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1]	434
[1320]	435	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1, &
	436	nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1, &
[1682]	437	nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) :: r !<
[1]	438
	439	!
	440	!-- Interpolate the residual
	441	l = grid_level
	442
[114]	443	!
	444	!-- Choose flag array of the upper level
[181]	445	SELECT CASE ( l+1 )
[114]	446	CASE ( 1 )
	447	flags => wall_flags_1
	448	CASE ( 2 )
	449	flags => wall_flags_2
	450	CASE ( 3 )
	451	flags => wall_flags_3
	452	CASE ( 4 )
	453	flags => wall_flags_4
	454	CASE ( 5 )
	455	flags => wall_flags_5
	456	CASE ( 6 )
	457	flags => wall_flags_6
	458	CASE ( 7 )
	459	flags => wall_flags_7
	460	CASE ( 8 )
	461	flags => wall_flags_8
	462	CASE ( 9 )
	463	flags => wall_flags_9
	464	CASE ( 10 )
	465	flags => wall_flags_10
	466	END SELECT
	467
[1]	468	!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
	469	!$OMP DO
	470	DO ic = nxl_mg(l), nxr_mg(l)
	471	i = 2*ic
	472	DO jc = nys_mg(l), nyn_mg(l)
	473	j = 2*jc
	474	DO kc = nzb+1, nzt_mg(l)
	475	k = 2*kc-1
[114]	476	!
	477	!-- Use implicit Neumann BCs if the respective gridpoint is inside
	478	!-- the building
	479	rkjim = r(k,j,i-1) + IBITS( flags(k,j,i-1), 6, 1 ) * &
	480	( r(k,j,i) - r(k,j,i-1) )
	481	rkjip = r(k,j,i+1) + IBITS( flags(k,j,i+1), 6, 1 ) * &
	482	( r(k,j,i) - r(k,j,i+1) )
	483	rkjpi = r(k,j+1,i) + IBITS( flags(k,j+1,i), 6, 1 ) * &
	484	( r(k,j,i) - r(k,j+1,i) )
	485	rkjmi = r(k,j-1,i) + IBITS( flags(k,j-1,i), 6, 1 ) * &
	486	( r(k,j,i) - r(k,j-1,i) )
	487	rkjmim = r(k,j-1,i-1) + IBITS( flags(k,j-1,i-1), 6, 1 ) * &
	488	( r(k,j,i) - r(k,j-1,i-1) )
	489	rkjpim = r(k,j+1,i-1) + IBITS( flags(k,j+1,i-1), 6, 1 ) * &
	490	( r(k,j,i) - r(k,j+1,i-1) )
	491	rkjmip = r(k,j-1,i+1) + IBITS( flags(k,j-1,i+1), 6, 1 ) * &
	492	( r(k,j,i) - r(k,j-1,i+1) )
	493	rkjpip = r(k,j+1,i+1) + IBITS( flags(k,j+1,i+1), 6, 1 ) * &
	494	( r(k,j,i) - r(k,j+1,i+1) )
	495	rkmji = r(k-1,j,i) + IBITS( flags(k-1,j,i), 6, 1 ) * &
	496	( r(k,j,i) - r(k-1,j,i) )
	497	rkmjim = r(k-1,j,i-1) + IBITS( flags(k-1,j,i-1), 6, 1 ) * &
	498	( r(k,j,i) - r(k-1,j,i-1) )
	499	rkmjip = r(k-1,j,i+1) + IBITS( flags(k-1,j,i+1), 6, 1 ) * &
	500	( r(k,j,i) - r(k-1,j,i+1) )
	501	rkmjpi = r(k-1,j+1,i) + IBITS( flags(k-1,j+1,i), 6, 1 ) * &
	502	( r(k,j,i) - r(k-1,j+1,i) )
	503	rkmjmi = r(k-1,j-1,i) + IBITS( flags(k-1,j-1,i), 6, 1 ) * &
	504	( r(k,j,i) - r(k-1,j-1,i) )
	505	rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
	506	( r(k,j,i) - r(k-1,j-1,i-1) )
	507	rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
	508	( r(k,j,i) - r(k-1,j+1,i-1) )
	509	rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
	510	( r(k,j,i) - r(k-1,j-1,i+1) )
	511	rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
	512	( r(k,j,i) - r(k-1,j+1,i+1) )
	513
[1353]	514	f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
	515	8.0_wp * r(k,j,i) &
	516	+ 4.0_wp * ( rkjim + rkjip + &
	517	rkjpi + rkjmi ) &
	518	+ 2.0_wp * ( rkjmim + rkjpim + &
	519	rkjmip + rkjpip ) &
	520	+ 4.0_wp * rkmji &
	521	+ 2.0_wp * ( rkmjim + rkmjim + &
	522	rkmjpi + rkmjmi ) &
	523	+ ( rkmjmim + rkmjpim + &
	524	rkmjmip + rkmjpip ) &
	525	+ 4.0_wp * r(k+1,j,i) &
	526	+ 2.0_wp * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
	527	r(k+1,j+1,i) + r(k+1,j-1,i) ) &
	528	+ ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
	529	r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
	530	)
[114]	531
[1353]	532	! f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * ( &
	533	! 8.0_wp * r(k,j,i) &
	534	! + 4.0_wp * ( r(k,j,i-1) + r(k,j,i+1) + &
	535	! r(k,j+1,i) + r(k,j-1,i) ) &
	536	! + 2.0_wp * ( r(k,j-1,i-1) + r(k,j+1,i-1) + &
	537	! r(k,j-1,i+1) + r(k,j+1,i+1) ) &
	538	! + 4.0_wp * r(k-1,j,i) &
	539	! + 2.0_wp * ( r(k-1,j,i-1) + r(k-1,j,i+1) + &
	540	! r(k-1,j+1,i) + r(k-1,j-1,i) ) &
	541	! + ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
	542	! r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
	543	! + 4.0_wp * r(k+1,j,i) &
	544	! + 2.0_wp * ( r(k+1,j,i-1) + r(k+1,j,i+1) + &
	545	! r(k+1,j+1,i) + r(k+1,j-1,i) ) &
	546	! + ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
	547	! r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
	548	! )
[1]	549	ENDDO
	550	ENDDO
	551	ENDDO
	552	!$OMP END PARALLEL
	553
	554	!
	555	!-- Horizontal boundary conditions
[667]	556	CALL exchange_horiz( f_mg, 1)
[1]	557
[707]	558	IF ( .NOT. bc_lr_cyc ) THEN
[1]	559	IF (inflow_l .OR. outflow_l) f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
	560	IF (inflow_r .OR. outflow_r) f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
	561	ENDIF
	562
[707]	563	IF ( .NOT. bc_ns_cyc ) THEN
[1]	564	IF (inflow_n .OR. outflow_n) f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
	565	IF (inflow_s .OR. outflow_s) f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
	566	ENDIF
	567
	568	!
[707]	569	!-- Boundary conditions at bottom and top of the domain.
	570	!-- These points are not handled by the above loop. Points may be within
	571	!-- buildings, but that doesn't matter.
	572	IF ( ibc_p_b == 1 ) THEN
	573	f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
	574	ELSE
[1353]	575	f_mg(nzb,:,: ) = 0.0_wp
[707]	576	ENDIF
[1]	577
[707]	578	IF ( ibc_p_t == 1 ) THEN
	579	f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
	580	ELSE
[1353]	581	f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
[707]	582	ENDIF
[1]	583
[707]	584
[1]	585	END SUBROUTINE restrict
	586
	587
	588	!------------------------------------------------------------------------------!
	589	! Description:
	590	! ------------
[1682]	591	!> Interpolates the correction of the perturbation pressure
	592	!> to the next finer grid.
[1]	593	!------------------------------------------------------------------------------!
[1682]	594	SUBROUTINE prolong( p, temp )
[1]	595
[1682]	596
[1320]	597	USE control_parameters, &
	598	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
	599	inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r, &
	600	outflow_s
[1]	601
[1320]	602	USE indices, &
	603	ONLY: nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	604
	605	USE kinds
	606
[1]	607	IMPLICIT NONE
	608
[1682]	609	INTEGER(iwp) :: i !<
	610	INTEGER(iwp) :: j !<
	611	INTEGER(iwp) :: k !<
	612	INTEGER(iwp) :: l !<
[1]	613
[1320]	614	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	615	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	616	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) :: p !<
[1]	617
[1320]	618	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	619	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	620	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: temp !<
[1]	621
	622
	623	!
	624	!-- First, store elements of the coarser grid on the next finer grid
	625	l = grid_level
	626
	627	!$OMP PARALLEL PRIVATE (i,j,k)
	628	!$OMP DO
	629	DO i = nxl_mg(l-1), nxr_mg(l-1)
	630	DO j = nys_mg(l-1), nyn_mg(l-1)
	631	!CDIR NODEP
	632	DO k = nzb+1, nzt_mg(l-1)
	633	!
	634	!-- Points of the coarse grid are directly stored on the next finer
	635	!-- grid
	636	temp(2k-1,2j,2*i) = p(k,j,i)
	637	!
	638	!-- Points between two coarse-grid points
[1353]	639	temp(2k-1,2j,2i+1) = 0.5_wp ( p(k,j,i) + p(k,j,i+1) )
	640	temp(2k-1,2j+1,2i) = 0.5_wp ( p(k,j,i) + p(k,j+1,i) )
	641	temp(2k,2j,2i) = 0.5_wp ( p(k,j,i) + p(k+1,j,i) )
[1]	642	!
	643	!-- Points in the center of the planes stretched by four points
	644	!-- of the coarse grid cube
[1353]	645	temp(2k-1,2j+1,2i+1) = 0.25_wp ( p(k,j,i) + p(k,j,i+1) + &
	646	p(k,j+1,i) + p(k,j+1,i+1) )
	647	temp(2k,2j,2i+1) = 0.25_wp ( p(k,j,i) + p(k,j,i+1) + &
	648	p(k+1,j,i) + p(k+1,j,i+1) )
	649	temp(2k,2j+1,2i) = 0.25_wp ( p(k,j,i) + p(k,j+1,i) + &
	650	p(k+1,j,i) + p(k+1,j+1,i) )
[1]	651	!
	652	!-- Points in the middle of coarse grid cube
[1353]	653	temp(2k,2j+1,2i+1) = 0.125_wp ( p(k,j,i) + p(k,j,i+1) + &
	654	p(k,j+1,i) + p(k,j+1,i+1) + &
	655	p(k+1,j,i) + p(k+1,j,i+1) + &
	656	p(k+1,j+1,i) + p(k+1,j+1,i+1) )
[1]	657	ENDDO
	658	ENDDO
	659	ENDDO
	660	!$OMP END PARALLEL
	661
	662	!
	663	!-- Horizontal boundary conditions
[667]	664	CALL exchange_horiz( temp, 1)
[1]	665
[707]	666	IF ( .NOT. bc_lr_cyc ) THEN
[1]	667	IF (inflow_l .OR. outflow_l) temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
	668	IF (inflow_r .OR. outflow_r) temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
	669	ENDIF
	670
[707]	671	IF ( .NOT. bc_ns_cyc ) THEN
[1]	672	IF (inflow_n .OR. outflow_n) temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
	673	IF (inflow_s .OR. outflow_s) temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
	674	ENDIF
	675
	676	!
	677	!-- Bottom and top boundary conditions
	678	IF ( ibc_p_b == 1 ) THEN
	679	temp(nzb,:,: ) = temp(nzb+1,:,:)
	680	ELSE
[1353]	681	temp(nzb,:,: ) = 0.0_wp
[1]	682	ENDIF
	683
	684	IF ( ibc_p_t == 1 ) THEN
	685	temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
	686	ELSE
[1353]	687	temp(nzt_mg(l)+1,:,: ) = 0.0_wp
[1]	688	ENDIF
	689
	690
	691	END SUBROUTINE prolong
	692
	693
	694	!------------------------------------------------------------------------------!
	695	! Description:
	696	! ------------
[1682]	697	!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
	698	!> 3D-Red-Black decomposition (GS-RB) is used.
[1]	699	!------------------------------------------------------------------------------!
[1682]	700	SUBROUTINE redblack( f_mg, p_mg )
[1]	701
[1682]	702
[1320]	703	USE arrays_3d, &
	704	ONLY: f1_mg, f2_mg, f3_mg
[1]	705
[1320]	706	USE control_parameters, &
	707	ONLY: bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l, &
	708	inflow_n, inflow_r, inflow_s, ngsrb, outflow_l, outflow_n, &
	709	outflow_r, outflow_s
	710
	711	USE cpulog, &
	712	ONLY: cpu_log, log_point_s
	713
	714	USE grid_variables, &
	715	ONLY: ddx2_mg, ddy2_mg
	716
	717	USE indices, &
	718	ONLY: flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4, &
	719	wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8, &
	720	wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg, &
	721	nzb, nzt_mg
	722
	723	USE kinds
	724
[1]	725	IMPLICIT NONE
	726
[1682]	727	INTEGER(iwp) :: color !<
	728	INTEGER(iwp) :: i !<
	729	INTEGER(iwp) :: ic !<
	730	INTEGER(iwp) :: j !<
	731	INTEGER(iwp) :: jc !<
	732	INTEGER(iwp) :: jj !<
	733	INTEGER(iwp) :: k !<
	734	INTEGER(iwp) :: l !<
	735	INTEGER(iwp) :: n !<
[1]	736
[1682]	737	LOGICAL :: unroll !<
[1]	738
[1682]	739	REAL(wp) :: wall_left !<
	740	REAL(wp) :: wall_north !<
	741	REAL(wp) :: wall_right !<
	742	REAL(wp) :: wall_south !<
	743	REAL(wp) :: wall_total !<
	744	REAL(wp) :: wall_top !<
[114]	745
[1320]	746	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	747	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	748	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1320]	749	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	750	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	751	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1]	752
	753	l = grid_level
	754
[114]	755	!
	756	!-- Choose flag array of this level
	757	SELECT CASE ( l )
	758	CASE ( 1 )
	759	flags => wall_flags_1
	760	CASE ( 2 )
	761	flags => wall_flags_2
	762	CASE ( 3 )
	763	flags => wall_flags_3
	764	CASE ( 4 )
	765	flags => wall_flags_4
	766	CASE ( 5 )
	767	flags => wall_flags_5
	768	CASE ( 6 )
	769	flags => wall_flags_6
	770	CASE ( 7 )
	771	flags => wall_flags_7
	772	CASE ( 8 )
	773	flags => wall_flags_8
	774	CASE ( 9 )
	775	flags => wall_flags_9
	776	CASE ( 10 )
	777	flags => wall_flags_10
	778	END SELECT
	779
[1]	780	unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0 .AND. &
	781	MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
	782
	783	DO n = 1, ngsrb
	784
[1320]	785	DO color = 1, 2
[1]	786
	787	IF ( .NOT. unroll ) THEN
[778]	788
[1]	789	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
	790
	791	!
	792	!-- Without unrolling of loops, no cache optimization
	793	DO i = nxl_mg(l), nxr_mg(l), 2
[1320]	794	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
[1]	795	DO k = nzb+1, nzt_mg(l), 2
[1353]	796	! p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	797	! ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
	798	! + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
	799	! + f2_mg(k,l) * p_mg(k+1,j,i) &
	800	! + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i) &
	801	! )
	802
[1353]	803	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	804	ddx2_mg(l) * &
	805	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	806	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	807	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	808	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	809	+ ddy2_mg(l) * &
	810	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	811	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	812	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	813	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	814	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	815	+ f3_mg(k,l) * &
	816	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	817	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	818	- f_mg(k,j,i) )
[1]	819	ENDDO
	820	ENDDO
	821	ENDDO
	822
	823	DO i = nxl_mg(l)+1, nxr_mg(l), 2
[1320]	824	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
[1]	825	DO k = nzb+1, nzt_mg(l), 2
[1353]	826	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	827	ddx2_mg(l) * &
	828	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	829	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	830	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	831	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	832	+ ddy2_mg(l) * &
	833	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	834	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	835	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	836	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	837	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	838	+ f3_mg(k,l) * &
	839	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	840	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	841	- f_mg(k,j,i) )
[1]	842	ENDDO
	843	ENDDO
	844	ENDDO
	845
	846	DO i = nxl_mg(l), nxr_mg(l), 2
[1320]	847	DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
[1]	848	DO k = nzb+2, nzt_mg(l), 2
[1353]	849	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	850	ddx2_mg(l) * &
	851	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	852	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	853	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	854	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	855	+ ddy2_mg(l) * &
	856	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	857	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	858	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	859	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	860	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	861	+ f3_mg(k,l) * &
	862	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	863	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	864	- f_mg(k,j,i) )
[1]	865	ENDDO
	866	ENDDO
	867	ENDDO
	868
	869	DO i = nxl_mg(l)+1, nxr_mg(l), 2
[1320]	870	DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
[1]	871	DO k = nzb+2, nzt_mg(l), 2
[1353]	872	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	873	ddx2_mg(l) * &
	874	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	875	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	876	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	877	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	878	+ ddy2_mg(l) * &
	879	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	880	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	881	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	882	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	883	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	884	+ f3_mg(k,l) * &
	885	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	886	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	887	- f_mg(k,j,i) )
[1]	888	ENDDO
	889	ENDDO
	890	ENDDO
	891	CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
	892
	893	ELSE
	894
	895	!
	896	!-- Loop unrolling along y, only one i loop for better cache use
	897	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
	898	DO ic = nxl_mg(l), nxr_mg(l), 2
	899	DO jc = nys_mg(l), nyn_mg(l), 4
	900	i = ic
[1320]	901	jj = jc+2-color
[1]	902	DO k = nzb+1, nzt_mg(l), 2
	903	j = jj
[1353]	904	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	905	ddx2_mg(l) * &
	906	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	907	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	908	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	909	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	910	+ ddy2_mg(l) * &
	911	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	912	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	913	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	914	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	915	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	916	+ f3_mg(k,l) * &
	917	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	918	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	919	- f_mg(k,j,i) )
[1]	920	j = jj+2
[1353]	921	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	922	ddx2_mg(l) * &
	923	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	924	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	925	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	926	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	927	+ ddy2_mg(l) * &
	928	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	929	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	930	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	931	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	932	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	933	+ f3_mg(k,l) * &
	934	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	935	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	936	- f_mg(k,j,i) )
[1]	937	ENDDO
	938
	939	i = ic+1
[1320]	940	jj = jc+color-1
[1]	941	DO k = nzb+1, nzt_mg(l), 2
	942	j =jj
[1353]	943	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	944	ddx2_mg(l) * &
	945	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	946	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	947	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	948	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	949	+ ddy2_mg(l) * &
	950	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	951	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	952	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	953	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	954	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	955	+ f3_mg(k,l) * &
	956	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	957	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	958	- f_mg(k,j,i) )
[1]	959	j = jj+2
[1353]	960	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	961	ddx2_mg(l) * &
	962	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	963	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	964	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	965	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	966	+ ddy2_mg(l) * &
	967	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	968	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	969	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	970	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	971	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	972	+ f3_mg(k,l) * &
	973	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	974	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	975	- f_mg(k,j,i) )
[1]	976	ENDDO
	977
	978	i = ic
[1320]	979	jj = jc+color-1
[1]	980	DO k = nzb+2, nzt_mg(l), 2
	981	j =jj
[1353]	982	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	983	ddx2_mg(l) * &
	984	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	985	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	986	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	987	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	988	+ ddy2_mg(l) * &
	989	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	990	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	991	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	992	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	993	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	994	+ f3_mg(k,l) * &
	995	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	996	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	997	- f_mg(k,j,i) )
[1]	998	j = jj+2
[1353]	999	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	1000	ddx2_mg(l) * &
	1001	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1002	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1003	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1004	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1005	+ ddy2_mg(l) * &
	1006	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1007	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1008	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1009	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1010	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	1011	+ f3_mg(k,l) * &
	1012	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1013	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1014	- f_mg(k,j,i) )
[1]	1015	ENDDO
	1016
	1017	i = ic+1
[1320]	1018	jj = jc+2-color
[1]	1019	DO k = nzb+2, nzt_mg(l), 2
	1020	j =jj
[1353]	1021	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	1022	ddx2_mg(l) * &
	1023	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1024	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1025	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1026	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1027	+ ddy2_mg(l) * &
	1028	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1029	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1030	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1031	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1032	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	1033	+ f3_mg(k,l) * &
	1034	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1035	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1036	- f_mg(k,j,i) )
[1]	1037	j = jj+2
[1353]	1038	p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * ( &
[114]	1039	ddx2_mg(l) * &
	1040	( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
	1041	( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
	1042	p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
	1043	( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
	1044	+ ddy2_mg(l) * &
	1045	( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
	1046	( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
	1047	p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
	1048	( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
	1049	+ f2_mg(k,l) * p_mg(k+1,j,i) &
	1050	+ f3_mg(k,l) * &
	1051	( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
	1052	( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
	1053	- f_mg(k,j,i) )
[1]	1054	ENDDO
	1055
	1056	ENDDO
	1057	ENDDO
	1058	CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
	1059
	1060	ENDIF
	1061
	1062	!
	1063	!-- Horizontal boundary conditions
[667]	1064	CALL exchange_horiz( p_mg, 1 )
[1]	1065
[707]	1066	IF ( .NOT. bc_lr_cyc ) THEN
[1]	1067	IF ( inflow_l .OR. outflow_l ) THEN
	1068	p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
	1069	ENDIF
	1070	IF ( inflow_r .OR. outflow_r ) THEN
	1071	p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
	1072	ENDIF
	1073	ENDIF
	1074
[707]	1075	IF ( .NOT. bc_ns_cyc ) THEN
[1]	1076	IF ( inflow_n .OR. outflow_n ) THEN
	1077	p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
	1078	ENDIF
	1079	IF ( inflow_s .OR. outflow_s ) THEN
	1080	p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
	1081	ENDIF
	1082	ENDIF
	1083
	1084	!
	1085	!-- Bottom and top boundary conditions
	1086	IF ( ibc_p_b == 1 ) THEN
	1087	p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
	1088	ELSE
[1353]	1089	p_mg(nzb,:,: ) = 0.0_wp
[1]	1090	ENDIF
	1091
	1092	IF ( ibc_p_t == 1 ) THEN
	1093	p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
	1094	ELSE
[1353]	1095	p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
[1]	1096	ENDIF
	1097
	1098	ENDDO
	1099
	1100	ENDDO
	1101
[114]	1102	!
	1103	!-- Set pressure within topography and at the topography surfaces
	1104	!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
	1105	!$OMP DO
	1106	DO i = nxl_mg(l), nxr_mg(l)
	1107	DO j = nys_mg(l), nyn_mg(l)
	1108	DO k = nzb, nzt_mg(l)
	1109	!
	1110	!-- First, set pressure inside topography to zero
[1353]	1111	p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
[114]	1112	!
	1113	!-- Second, determine if the gridpoint inside topography is adjacent
	1114	!-- to a wall and set its value to a value given by the average of
	1115	!-- those values obtained from Neumann boundary condition
	1116	wall_left = IBITS( flags(k,j,i-1), 5, 1 )
	1117	wall_right = IBITS( flags(k,j,i+1), 4, 1 )
	1118	wall_south = IBITS( flags(k,j-1,i), 3, 1 )
	1119	wall_north = IBITS( flags(k,j+1,i), 2, 1 )
	1120	wall_top = IBITS( flags(k+1,j,i), 0, 1 )
	1121	wall_total = wall_left + wall_right + wall_south + wall_north + &
	1122	wall_top
[1]	1123
[1353]	1124	IF ( wall_total > 0.0_wp ) THEN
	1125	p_mg(k,j,i) = 1.0_wp / wall_total * &
	1126	( wall_left * p_mg(k,j,i-1) + &
	1127	wall_right * p_mg(k,j,i+1) + &
	1128	wall_south * p_mg(k,j-1,i) + &
	1129	wall_north * p_mg(k,j+1,i) + &
	1130	wall_top * p_mg(k+1,j,i) )
[114]	1131	ENDIF
	1132	ENDDO
	1133	ENDDO
	1134	ENDDO
[1056]	1135	!$OMP END PARALLEL
[114]	1136
	1137	!
	1138	!-- One more time horizontal boundary conditions
[667]	1139	CALL exchange_horiz( p_mg, 1)
[114]	1140
[778]	1141
[1]	1142	END SUBROUTINE redblack
	1143
	1144
	1145
[1682]	1146	!------------------------------------------------------------------------------!
	1147	! Description:
	1148	! ------------
	1149	!> Gather subdomain data from all PEs.
	1150	!------------------------------------------------------------------------------!
[1]	1151	SUBROUTINE mg_gather( f2, f2_sub )
	1152
[1320]	1153	USE control_parameters, &
	1154	ONLY: grid_level
	1155
	1156	USE cpulog, &
	1157	ONLY: cpu_log, log_point_s
	1158
	1159	USE indices, &
	1160	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1161
	1162	USE kinds
	1163
[1]	1164	USE pegrid
	1165
	1166	IMPLICIT NONE
	1167
[1682]	1168	INTEGER(iwp) :: i !<
	1169	INTEGER(iwp) :: il !<
	1170	INTEGER(iwp) :: ir !<
	1171	INTEGER(iwp) :: j !<
	1172	INTEGER(iwp) :: jn !<
	1173	INTEGER(iwp) :: js !<
	1174	INTEGER(iwp) :: k !<
	1175	INTEGER(iwp) :: nwords !<
[1]	1176
[1320]	1177	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1178	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1179	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2 !<
[1320]	1180	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1181	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1182	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f2_l !<
[1]	1183
[1320]	1184	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
	1185	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
[1682]	1186	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: f2_sub !<
[1]	1187
	1188
	1189	#if defined( __parallel )
	1190	CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
	1191
[1353]	1192	f2_l = 0.0_wp
[1]	1193
	1194	!
[707]	1195	!-- Store the local subdomain array on the total array
	1196	js = mg_loc_ind(3,myid)
	1197	IF ( south_border_pe ) js = js - 1
	1198	jn = mg_loc_ind(4,myid)
	1199	IF ( north_border_pe ) jn = jn + 1
	1200	il = mg_loc_ind(1,myid)
	1201	IF ( left_border_pe ) il = il - 1
	1202	ir = mg_loc_ind(2,myid)
	1203	IF ( right_border_pe ) ir = ir + 1
	1204	DO i = il, ir
	1205	DO j = js, jn
	1206	DO k = nzb, nzt_mg(grid_level)+1
	1207	f2_l(k,j,i) = f2_sub(k,j,i)
	1208	ENDDO
[1]	1209	ENDDO
[707]	1210	ENDDO
[1]	1211
	1212	!
[707]	1213	!-- Find out the number of array elements of the total array
	1214	nwords = SIZE( f2 )
[1]	1215
[707]	1216	!
	1217	!-- Gather subdomain data from all PEs
	1218	IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
	1219	CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
	1220	f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
	1221	nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
	1222
[1]	1223	CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
	1224	#endif
	1225
	1226	END SUBROUTINE mg_gather
	1227
	1228
	1229
[1682]	1230	!------------------------------------------------------------------------------!
	1231	! Description:
	1232	! ------------
	1233	!> @todo It may be possible to improve the speed of this routine by using
	1234	!> non-blocking communication
	1235	!------------------------------------------------------------------------------!
[1]	1236	SUBROUTINE mg_scatter( p2, p2_sub )
	1237
[1320]	1238	USE control_parameters, &
	1239	ONLY: grid_level
	1240
	1241	USE cpulog, &
	1242	ONLY: cpu_log, log_point_s
	1243
	1244	USE indices, &
	1245	ONLY: mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
	1246
	1247	USE kinds
	1248
[1]	1249	USE pegrid
	1250
	1251	IMPLICIT NONE
	1252
[1682]	1253	INTEGER(iwp) :: nwords !<
[1]	1254
[1320]	1255	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1256	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1257	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[1]	1258
[1320]	1259	REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1, &
	1260	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
[1682]	1261	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) :: p2_sub !<
[1]	1262
	1263	!
	1264	!-- Find out the number of array elements of the subdomain array
	1265	nwords = SIZE( p2_sub )
	1266
	1267	#if defined( __parallel )
	1268	CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
	1269
[707]	1270	p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1271	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
[1]	1272
	1273	CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
	1274	#endif
	1275
	1276	END SUBROUTINE mg_scatter
	1277
	1278
	1279	!------------------------------------------------------------------------------!
	1280	! Description:
	1281	! ------------
[1682]	1282	!> This is where the multigrid technique takes place. V- and W- Cycle are
	1283	!> implemented and steered by the parameter "gamma". Parameter "nue" determines
	1284	!> the convergence of the multigrid iterative solution. There are nue times
	1285	!> RB-GS iterations. It should be set to "1" or "2", considering the time effort
	1286	!> one would like to invest. Last choice shows a very good converging factor,
	1287	!> but leads to an increase in computing time.
[1]	1288	!------------------------------------------------------------------------------!
[1682]	1289	RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
[1]	1290
[1320]	1291	USE control_parameters, &
	1292	ONLY: bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir, &
	1293	gamma_mg, grid_level, grid_level_count, ibc_p_b, ibc_p_t, &
	1294	inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level, &
	1295	mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb, outflow_l, &
	1296	outflow_n, outflow_r, outflow_s
	1297
	1298
	1299	USE indices, &
	1300	ONLY: mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn, &
	1301	nyn_mg, nzb, nzt, nzt_mg
	1302
	1303	USE kinds
	1304
[1]	1305	USE pegrid
	1306
	1307	IMPLICIT NONE
	1308
[1682]	1309	INTEGER(iwp) :: i !<
	1310	INTEGER(iwp) :: j !<
	1311	INTEGER(iwp) :: k !<
	1312	INTEGER(iwp) :: nxl_mg_save !<
	1313	INTEGER(iwp) :: nxr_mg_save !<
	1314	INTEGER(iwp) :: nyn_mg_save !<
	1315	INTEGER(iwp) :: nys_mg_save !<
	1316	INTEGER(iwp) :: nzt_mg_save !<
[1]	1317
[1320]	1318	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1319	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1320	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg !<
[1320]	1321	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1322	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1323	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg !<
[1320]	1324	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1325	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1326	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3 !<
[1320]	1327	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1, &
	1328	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
[1682]	1329	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r !<
[1]	1330
[1320]	1331	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1332	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1333	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: f2 !<
[1320]	1334	REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1, &
	1335	nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1, &
[1682]	1336	nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) :: p2 !<
[1]	1337
[1682]	1338	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: f2_sub !<
	1339	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: p2_sub !<
[778]	1340
[1]	1341	!
	1342	!-- Restriction to the coarsest grid
	1343	10 IF ( grid_level == 1 ) THEN
	1344
	1345	!
	1346	!-- Solution on the coarsest grid. Double the number of Gauss-Seidel
	1347	!-- iterations in order to get a more accurate solution.
	1348	ngsrb = 2 * ngsrb
[778]	1349
[1]	1350	CALL redblack( f_mg, p_mg )
[778]	1351
[1]	1352	ngsrb = ngsrb / 2
	1353
[778]	1354
[1]	1355	ELSEIF ( grid_level /= 1 ) THEN
	1356
	1357	grid_level_count(grid_level) = grid_level_count(grid_level) + 1
	1358
	1359	!
	1360	!-- Solution on the actual grid level
	1361	CALL redblack( f_mg, p_mg )
	1362
	1363	!
	1364	!-- Determination of the actual residual
	1365	CALL resid( f_mg, p_mg, r )
	1366
	1367	!
	1368	!-- Restriction of the residual (finer grid values!) to the next coarser
	1369	!-- grid. Therefore, the grid level has to be decremented now. nxl..nzt have
	1370	!-- to be set to the coarse grid values, because these variables are needed
	1371	!-- for the exchange of ghost points in routine exchange_horiz
	1372	grid_level = grid_level - 1
	1373	nxl = nxl_mg(grid_level)
[778]	1374	nys = nys_mg(grid_level)
[1]	1375	nxr = nxr_mg(grid_level)
	1376	nyn = nyn_mg(grid_level)
	1377	nzt = nzt_mg(grid_level)
	1378
	1379	IF ( grid_level == mg_switch_to_pe0_level ) THEN
[778]	1380
[1]	1381	!
	1382	!-- From this level on, calculations are done on PE0 only.
	1383	!-- First, carry out restriction on the subdomain.
	1384	!-- Therefore, indices of the level have to be changed to subdomain values
	1385	!-- in between (otherwise, the restrict routine would expect
	1386	!-- the gathered array)
[778]	1387
[1]	1388	nxl_mg_save = nxl_mg(grid_level)
	1389	nxr_mg_save = nxr_mg(grid_level)
	1390	nys_mg_save = nys_mg(grid_level)
	1391	nyn_mg_save = nyn_mg(grid_level)
	1392	nzt_mg_save = nzt_mg(grid_level)
	1393	nxl_mg(grid_level) = mg_loc_ind(1,myid)
	1394	nxr_mg(grid_level) = mg_loc_ind(2,myid)
	1395	nys_mg(grid_level) = mg_loc_ind(3,myid)
	1396	nyn_mg(grid_level) = mg_loc_ind(4,myid)
	1397	nzt_mg(grid_level) = mg_loc_ind(5,myid)
	1398	nxl = mg_loc_ind(1,myid)
	1399	nxr = mg_loc_ind(2,myid)
	1400	nys = mg_loc_ind(3,myid)
	1401	nyn = mg_loc_ind(4,myid)
	1402	nzt = mg_loc_ind(5,myid)
	1403
	1404	ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1, &
	1405	nys_mg(grid_level)-1:nyn_mg(grid_level)+1, &
	1406	nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
	1407
	1408	CALL restrict( f2_sub, r )
	1409
	1410	!
	1411	!-- Restore the correct indices of this level
	1412	nxl_mg(grid_level) = nxl_mg_save
	1413	nxr_mg(grid_level) = nxr_mg_save
	1414	nys_mg(grid_level) = nys_mg_save
	1415	nyn_mg(grid_level) = nyn_mg_save
	1416	nzt_mg(grid_level) = nzt_mg_save
	1417	nxl = nxl_mg(grid_level)
	1418	nxr = nxr_mg(grid_level)
	1419	nys = nys_mg(grid_level)
	1420	nyn = nyn_mg(grid_level)
	1421	nzt = nzt_mg(grid_level)
	1422	!
	1423	!-- Gather all arrays from the subdomains on PE0
	1424	CALL mg_gather( f2, f2_sub )
	1425
	1426	!
	1427	!-- Set switch for routine exchange_horiz, that no ghostpoint exchange
	1428	!-- has to be carried out from now on
	1429	mg_switch_to_pe0 = .TRUE.
	1430
	1431	!
	1432	!-- In case of non-cyclic lateral boundary conditions, both in- and
[707]	1433	!-- outflow conditions have to be used on all PEs after the switch,
	1434	!-- because then they have the total domain.
[1159]	1435	IF ( bc_lr_dirrad ) THEN
[707]	1436	inflow_l = .TRUE.
	1437	inflow_r = .FALSE.
	1438	outflow_l = .FALSE.
	1439	outflow_r = .TRUE.
[1159]	1440	ELSEIF ( bc_lr_raddir ) THEN
[707]	1441	inflow_l = .FALSE.
	1442	inflow_r = .TRUE.
	1443	outflow_l = .TRUE.
	1444	outflow_r = .FALSE.
[1]	1445	ENDIF
	1446
[1159]	1447	IF ( bc_ns_dirrad ) THEN
[707]	1448	inflow_n = .TRUE.
	1449	inflow_s = .FALSE.
	1450	outflow_n = .FALSE.
	1451	outflow_s = .TRUE.
[1159]	1452	ELSEIF ( bc_ns_raddir ) THEN
[707]	1453	inflow_n = .FALSE.
	1454	inflow_s = .TRUE.
	1455	outflow_n = .TRUE.
	1456	outflow_s = .FALSE.
	1457	ENDIF
	1458
[1]	1459	DEALLOCATE( f2_sub )
	1460
	1461	ELSE
[1056]	1462
[1]	1463	CALL restrict( f2, r )
	1464
	1465	ENDIF
[707]	1466
[1353]	1467	p2 = 0.0_wp
[1]	1468
	1469	!
	1470	!-- Repeat the same procedure till the coarsest grid is reached
[707]	1471	CALL next_mg_level( f2, p2, p3, r )
[1]	1472
	1473	ENDIF
	1474
	1475	!
	1476	!-- Now follows the prolongation
	1477	IF ( grid_level >= 2 ) THEN
	1478
	1479	!
	1480	!-- Prolongation of the new residual. The values are transferred
	1481	!-- from the coarse to the next finer grid.
	1482	IF ( grid_level == mg_switch_to_pe0_level+1 ) THEN
[879]	1483
	1484	#if defined( __parallel )
[1]	1485	!
	1486	!-- At this level, the new residual first has to be scattered from
	1487	!-- PE0 to the other PEs
	1488	ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1, &
	1489	mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
	1490	mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
	1491
	1492	CALL mg_scatter( p2, p2_sub )
	1493
	1494	!
	1495	!-- Therefore, indices of the previous level have to be changed to
	1496	!-- subdomain values in between (otherwise, the prolong routine would
	1497	!-- expect the gathered array)
	1498	nxl_mg_save = nxl_mg(grid_level-1)
	1499	nxr_mg_save = nxr_mg(grid_level-1)
	1500	nys_mg_save = nys_mg(grid_level-1)
	1501	nyn_mg_save = nyn_mg(grid_level-1)
	1502	nzt_mg_save = nzt_mg(grid_level-1)
	1503	nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
	1504	nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
	1505	nys_mg(grid_level-1) = mg_loc_ind(3,myid)
	1506	nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
	1507	nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
	1508
	1509	!
	1510	!-- Set switch for routine exchange_horiz, that ghostpoint exchange
	1511	!-- has to be carried again out from now on
	1512	mg_switch_to_pe0 = .FALSE.
	1513
	1514	!
[707]	1515	!-- For non-cyclic lateral boundary conditions, restore the
	1516	!-- in-/outflow conditions
	1517	inflow_l = .FALSE.; inflow_r = .FALSE.
	1518	inflow_n = .FALSE.; inflow_s = .FALSE.
	1519	outflow_l = .FALSE.; outflow_r = .FALSE.
	1520	outflow_n = .FALSE.; outflow_s = .FALSE.
	1521
	1522	IF ( pleft == MPI_PROC_NULL ) THEN
[1159]	1523	IF ( bc_lr_dirrad ) THEN
[707]	1524	inflow_l = .TRUE.
[1159]	1525	ELSEIF ( bc_lr_raddir ) THEN
[707]	1526	outflow_l = .TRUE.
[1]	1527	ENDIF
[707]	1528	ENDIF
	1529
	1530	IF ( pright == MPI_PROC_NULL ) THEN
[1159]	1531	IF ( bc_lr_dirrad ) THEN
[707]	1532	outflow_r = .TRUE.
[1159]	1533	ELSEIF ( bc_lr_raddir ) THEN
[707]	1534	inflow_r = .TRUE.
[1]	1535	ENDIF
	1536	ENDIF
	1537
[707]	1538	IF ( psouth == MPI_PROC_NULL ) THEN
[1159]	1539	IF ( bc_ns_dirrad ) THEN
[707]	1540	outflow_s = .TRUE.
[1159]	1541	ELSEIF ( bc_ns_raddir ) THEN
[707]	1542	inflow_s = .TRUE.
	1543	ENDIF
	1544	ENDIF
	1545
	1546	IF ( pnorth == MPI_PROC_NULL ) THEN
[1159]	1547	IF ( bc_ns_dirrad ) THEN
[707]	1548	inflow_n = .TRUE.
[1159]	1549	ELSEIF ( bc_ns_raddir ) THEN
[707]	1550	outflow_n = .TRUE.
	1551	ENDIF
	1552	ENDIF
	1553
[1]	1554	CALL prolong( p2_sub, p3 )
	1555
	1556	!
	1557	!-- Restore the correct indices of the previous level
	1558	nxl_mg(grid_level-1) = nxl_mg_save
	1559	nxr_mg(grid_level-1) = nxr_mg_save
	1560	nys_mg(grid_level-1) = nys_mg_save
	1561	nyn_mg(grid_level-1) = nyn_mg_save
	1562	nzt_mg(grid_level-1) = nzt_mg_save
	1563
	1564	DEALLOCATE( p2_sub )
[879]	1565	#endif
[1]	1566
	1567	ELSE
[879]	1568
[1]	1569	CALL prolong( p2, p3 )
	1570
	1571	ENDIF
	1572
	1573	!
	1574	!-- Computation of the new pressure correction. Therefore,
	1575	!-- values from prior grids are added up automatically stage by stage.
	1576	DO i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
	1577	DO j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
	1578	DO k = nzb, nzt_mg(grid_level)+1
	1579	p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
	1580	ENDDO
	1581	ENDDO
	1582	ENDDO
	1583
	1584	!
	1585	!-- Relaxation of the new solution
	1586	CALL redblack( f_mg, p_mg )
	1587
	1588	ENDIF
	1589
[778]	1590
[1]	1591	!
	1592	!-- The following few lines serve the steering of the multigrid scheme
	1593	IF ( grid_level == maximum_grid_level ) THEN
	1594
	1595	GOTO 20
	1596
	1597	ELSEIF ( grid_level /= maximum_grid_level .AND. grid_level /= 1 .AND. &
	1598	grid_level_count(grid_level) /= gamma_mg ) THEN
	1599
	1600	GOTO 10
	1601
	1602	ENDIF
	1603
	1604	!
	1605	!-- Reset counter for the next call of poismg
	1606	grid_level_count(grid_level) = 0
	1607
	1608	!
	1609	!-- Continue with the next finer level. nxl..nzt have to be
	1610	!-- set to the finer grid values, because these variables are needed for the
	1611	!-- exchange of ghost points in routine exchange_horiz
	1612	grid_level = grid_level + 1
	1613	nxl = nxl_mg(grid_level)
	1614	nxr = nxr_mg(grid_level)
	1615	nys = nys_mg(grid_level)
	1616	nyn = nyn_mg(grid_level)
	1617	nzt = nzt_mg(grid_level)
	1618
	1619	20 CONTINUE
	1620
	1621	END SUBROUTINE next_mg_level

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |