Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |