Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt.f90 @ 2232

Last change on this file since 2232 was 2232, checked in by suehring, 8 years ago
Adjustments according new topography and surface-modelling concept implemented
Property svn:keywords set to `Id`
File size: 68.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |