Home

Context Navigation

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

Last change on this file since 251 was 198, checked in by raasch, 16 years ago
file headers updated for the next release 3.5
Property svn:keywords set to `Id`
File size: 54.2 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |