Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |