Home

Context Navigation

source: palm/trunk/SOURCE/poismg_noopt.f90 @ 2232

Last change on this file since 2232 was 2232, checked in by suehring, 7 years ago
Adjustments according new topography and surface-modelling concept implemented
Property svn:keywords set to `Id`
File size: 68.0 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |