Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |