Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |