Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |