Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |