Home

Context Navigation

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

Last change on this file since 883 was 881, checked in by raasch, 13 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 55.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |