Home

Context Navigation

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

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