Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |