Home

Context Navigation

poismg_noopt_mod.f90 @ 4844

Last change on this file since 4844 was 4832, checked in by raasch, 4 years ago
bugfix in redblack algorithm: lower i,j indices need to start alternatively with even or odd value on the coarsest grid level, if the subdomain has an uneven number of gridpoints along x/y; some debug output flushed
Property svn:keywords set to `Id`
File size: 96.6 KB

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

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

Context Navigation

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

Download in other formats: