Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |