Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |