Home

Context Navigation

source: palm/trunk/SOURCE/poismg.f90 @ 1267

Last change on this file since 1267 was 1160, checked in by fricke, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 56.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |