Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt_mod.f90 @ 2750

Last change on this file since 2750 was 2718, checked in by maronga, 7 years ago
deleting of deprecated files; headers updated where needed
Property svn:keywords set to `Id`
File size: 81.1 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |