Home

Context Navigation

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

Last change on this file since 1940 was 1935, checked in by hellstea, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 67.3 KB

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