Home

Context Navigation

poismg.f90 @ 707

Last change on this file since 707 was 707, checked in by raasch, 13 years ago

New:
---

In case of multigrid method, on coarse grid levels, gathered data are
identically processed on all PEs (before, on PE0 only), so that the subsequent
scattering of data is not neccessary any more. (modules, init_pegrid, poismg)

Changed:

Calculation of weighted average of p is now handled in the same way
regardless of the number of ghost layers (advection scheme). (pres)

multigrid and sor method are using p_loc for iterative
advancements of pressure. p_sub removed. (init_3d_model, modules, poismg, pres, sor)

bc_lr and bc_ns replaced by bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir
for speed optimization. (calc_spectra, check_parameters, exchange_horiz,
exchange_horiz_2d, header, init_3d_model, init_grid, init_pegrid, modules,
poismg, pres, sor, time_integration, timestep)

grid_level directly used as index for MPI data type arrays. (exchange_horiz,
poismg)

initial assignments of zero to array p for iterative solvers only (init_3d_model)

Errors:

localsum calculation modified for proper OpenMP reduction. (pres)

Bugfix: bottom (nzb) and top (nzt+1) boundary conditions set in routines
resid and restrict. They were missed before, which may have led to
unpredictable results. (poismg)

Property svn:keywords set to Id

File size: 54.6 KB

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

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

Context Navigation

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

Download in other formats: