Home

Context Navigation

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

Last change on this file since 697 was 668, checked in by suehring, 14 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 54.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |