Home

Context Navigation

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

Last change on this file since 4180 was 4180, checked in by scharf, 5 years ago
removed comments in 'Former revisions' section that are older than 01.01.2019
Property svn:keywords set to `Id`
File size: 76.8 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |