Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |