Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |