Home

Context Navigation

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

Last change on this file since 1041 was 1037, checked in by raasch, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 56.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |