Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |