Home

Context Navigation

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

Last change on this file since 1320 was 1320, checked in by raasch, 10 years ago
ONLY-attribute added to USE-statements, kind-parameters added to all INTEGER and REAL declaration statements, kinds are defined in new module kinds, old module precision_kind is removed, revision history before 2012 removed, comment fields (!:) to be used for variable explanations added to all variable declaration statements
Property svn:keywords set to `Id`
File size: 64.3 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |