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

Last change on this file since 1762 was 1762, checked in by hellstea, 6 years ago

Introduction of nested domain system

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