Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |