Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |