Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |