Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |