Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |