Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |