Home

Context Navigation

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

Last change on this file since 70 was 4, checked in by raasch, 18 years ago
Id keyword set as property for all *.f90 files
Property svn:keywords set to `Id`
File size: 39.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |