Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |