Home

Context Navigation

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

Last change on this file since 4628 was 4457, checked in by raasch, 5 years ago
ghost point exchange modularized, bugfix for wrong 2d-exchange
Property svn:keywords set to `Id`
File size: 77.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |