Home

Context Navigation

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

Last change on this file since 4235 was 4182, checked in by scharf, 5 years ago
corrected "Former revisions" section minor formatting in "Former revisions" section added "Author" section
Property svn:keywords set to `Id`
File size: 76.9 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |