Home

Context Navigation

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

Last change on this file since 4385 was 4360, checked in by suehring, 5 years ago
Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update
Property svn:keywords set to `Id`
File size: 77.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |