Home

Context Navigation

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

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