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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 65.3 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!
22!
23! Former revisions:
24! -----------------
25! $Id: poismg.f90 1683 2015-10-07 23:57:51Z knoop $
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, outflow_l, outflow_n, outflow_r,  &
262               outflow_s
263
264    USE grid_variables,                                                        &
265        ONLY:  ddx2_mg, ddy2_mg
266
267    USE indices,                                                               &
268        ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4,  &
269               wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,         &
270               wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg,    &
271               nzb, nzt_mg
272
273    USE kinds
274
275    IMPLICIT NONE
276
277    INTEGER(iwp) ::  i
278    INTEGER(iwp) ::  j
279    INTEGER(iwp) ::  k
280    INTEGER(iwp) ::  l
281
282    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
283                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
284                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
285    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
286                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
287                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
288    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
289                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
290                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !<
291
292!
293!-- Calculate the residual
294    l = grid_level
295
296!
297!-- Choose flag array of this level
298    SELECT CASE ( l )
299       CASE ( 1 )
300          flags => wall_flags_1
301       CASE ( 2 )
302          flags => wall_flags_2
303       CASE ( 3 )
304          flags => wall_flags_3
305       CASE ( 4 )
306          flags => wall_flags_4
307       CASE ( 5 )
308          flags => wall_flags_5
309       CASE ( 6 )
310          flags => wall_flags_6
311       CASE ( 7 )
312          flags => wall_flags_7
313       CASE ( 8 )
314          flags => wall_flags_8
315       CASE ( 9 )
316          flags => wall_flags_9
317       CASE ( 10 )
318          flags => wall_flags_10
319    END SELECT
320
321!$OMP PARALLEL PRIVATE (i,j,k)
322!$OMP DO
323    DO  i = nxl_mg(l), nxr_mg(l)
324       DO  j = nys_mg(l), nyn_mg(l) 
325          DO  k = nzb+1, nzt_mg(l)
326             r(k,j,i) = f_mg(k,j,i)                                         &
327                        - ddx2_mg(l) *                                      &
328                            ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
329                                          ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
330                              p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
331                                          ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
332                        - ddy2_mg(l) *                                      &
333                            ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
334                                          ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
335                              p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
336                                          ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
337                        - f2_mg(k,l) * p_mg(k+1,j,i)                        &
338                        - f3_mg(k,l) *                                      &
339                            ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
340                                          ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
341                        + f1_mg(k,l) * p_mg(k,j,i)
342!
343!--          Residual within topography should be zero
344             r(k,j,i) = r(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
345          ENDDO
346       ENDDO
347    ENDDO
348!$OMP END PARALLEL
349
350!
351!-- Horizontal boundary conditions
352    CALL exchange_horiz( r, 1)
353
354    IF ( .NOT. bc_lr_cyc )  THEN
355       IF ( inflow_l .OR. outflow_l )  r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
356       IF ( inflow_r .OR. outflow_r )  r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
357    ENDIF
358
359    IF ( .NOT. bc_ns_cyc )  THEN
360       IF ( inflow_n .OR. outflow_n )  r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
361       IF ( inflow_s .OR. outflow_s )  r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
362    ENDIF
363
364!
365!-- Boundary conditions at bottom and top of the domain.
366!-- These points are not handled by the above loop. Points may be within
367!-- buildings, but that doesn't matter.
368    IF ( ibc_p_b == 1 )  THEN
369       r(nzb,:,: ) = r(nzb+1,:,:)
370    ELSE
371       r(nzb,:,: ) = 0.0_wp
372    ENDIF
373
374    IF ( ibc_p_t == 1 )  THEN
375       r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
376    ELSE
377       r(nzt_mg(l)+1,:,: ) = 0.0_wp
378    ENDIF
379
380
381 END SUBROUTINE resid
382
383
384!------------------------------------------------------------------------------!
385! Description:
386! ------------
387!> Interpolates the residual on the next coarser grid with "full weighting"
388!> scheme.
389!------------------------------------------------------------------------------!
390 SUBROUTINE restrict( f_mg, r )
391
392
393    USE control_parameters,                                                    &
394        ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,   &
395               inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r,  &
396               outflow_s
397
398    USE indices,                                                               &
399        ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4,  &
400               wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,         &
401               wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg,    &
402               nzb, nzt_mg
403
404    USE kinds
405
406    IMPLICIT NONE
407
408    INTEGER(iwp) ::  i    !<
409    INTEGER(iwp) ::  ic   !<
410    INTEGER(iwp) ::  j    !<
411    INTEGER(iwp) ::  jc   !<
412    INTEGER(iwp) ::  k    !<
413    INTEGER(iwp) ::  kc   !<
414    INTEGER(iwp) ::  l    !<
415
416    REAL(wp) ::  rkjim    !<
417    REAL(wp) ::  rkjip    !<
418    REAL(wp) ::  rkjmi    !<
419    REAL(wp) ::  rkjmim   !<
420    REAL(wp) ::  rkjmip   !<
421    REAL(wp) ::  rkjpi    !<
422    REAL(wp) ::  rkjpim   !<
423    REAL(wp) ::  rkjpip   !<
424    REAL(wp) ::  rkmji    !<
425    REAL(wp) ::  rkmjim   !<
426    REAL(wp) ::  rkmjip   !<
427    REAL(wp) ::  rkmjmi   !<
428    REAL(wp) ::  rkmjmim  !<
429    REAL(wp) ::  rkmjmip  !<
430    REAL(wp) ::  rkmjpi   !<
431    REAL(wp) ::  rkmjpim  !<
432    REAL(wp) ::  rkmjpip  !<
433
434    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                            &
435                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,           &
436                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
437
438    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,                          &
439                        nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,       &
440                        nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r !<
441
442!
443!-- Interpolate the residual
444    l = grid_level
445
446!
447!-- Choose flag array of the upper level
448    SELECT CASE ( l+1 )
449       CASE ( 1 )
450          flags => wall_flags_1
451       CASE ( 2 )
452          flags => wall_flags_2
453       CASE ( 3 )
454          flags => wall_flags_3
455       CASE ( 4 )
456          flags => wall_flags_4
457       CASE ( 5 )
458          flags => wall_flags_5
459       CASE ( 6 )
460          flags => wall_flags_6
461       CASE ( 7 )
462          flags => wall_flags_7
463       CASE ( 8 )
464          flags => wall_flags_8
465       CASE ( 9 )
466          flags => wall_flags_9
467       CASE ( 10 )
468          flags => wall_flags_10
469    END SELECT
470   
471!$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc)
472!$OMP DO
473    DO  ic = nxl_mg(l), nxr_mg(l)   
474       i = 2*ic
475       DO  jc = nys_mg(l), nyn_mg(l) 
476          j = 2*jc
477          DO  kc = nzb+1, nzt_mg(l)
478             k = 2*kc-1
479!
480!--          Use implicit Neumann BCs if the respective gridpoint is inside
481!--          the building
482             rkjim   = r(k,j,i-1)     + IBITS( flags(k,j,i-1), 6, 1 ) *     &
483                                        ( r(k,j,i) - r(k,j,i-1) )
484             rkjip   = r(k,j,i+1)     + IBITS( flags(k,j,i+1), 6, 1 ) *     &
485                                        ( r(k,j,i) - r(k,j,i+1) )
486             rkjpi   = r(k,j+1,i)     + IBITS( flags(k,j+1,i), 6, 1 ) *     &
487                                        ( r(k,j,i) - r(k,j+1,i) )
488             rkjmi   = r(k,j-1,i)     + IBITS( flags(k,j-1,i), 6, 1 ) *     &
489                                        ( r(k,j,i) - r(k,j-1,i) )
490             rkjmim  = r(k,j-1,i-1)   + IBITS( flags(k,j-1,i-1), 6, 1 ) *   &
491                                        ( r(k,j,i) - r(k,j-1,i-1) )
492             rkjpim  = r(k,j+1,i-1)   + IBITS( flags(k,j+1,i-1), 6, 1 ) *   &
493                                        ( r(k,j,i) - r(k,j+1,i-1) )
494             rkjmip  = r(k,j-1,i+1)   + IBITS( flags(k,j-1,i+1), 6, 1 ) *   &
495                                        ( r(k,j,i) - r(k,j-1,i+1) )
496             rkjpip  = r(k,j+1,i+1)   + IBITS( flags(k,j+1,i+1), 6, 1 ) *   &
497                                        ( r(k,j,i) - r(k,j+1,i+1) )
498             rkmji   = r(k-1,j,i)     + IBITS( flags(k-1,j,i), 6, 1 ) *     &
499                                        ( r(k,j,i) - r(k-1,j,i) )
500             rkmjim  = r(k-1,j,i-1)   + IBITS( flags(k-1,j,i-1), 6, 1 ) *   &
501                                        ( r(k,j,i) - r(k-1,j,i-1) )
502             rkmjip  = r(k-1,j,i+1)   + IBITS( flags(k-1,j,i+1), 6, 1 ) *   &
503                                        ( r(k,j,i) - r(k-1,j,i+1) )
504             rkmjpi  = r(k-1,j+1,i)   + IBITS( flags(k-1,j+1,i), 6, 1 ) *   &
505                                        ( r(k,j,i) - r(k-1,j+1,i) )
506             rkmjmi  = r(k-1,j-1,i)   + IBITS( flags(k-1,j-1,i), 6, 1 ) *   &
507                                        ( r(k,j,i) - r(k-1,j-1,i) )
508             rkmjmim = r(k-1,j-1,i-1) + IBITS( flags(k-1,j-1,i-1), 6, 1 ) * &
509                                        ( r(k,j,i) - r(k-1,j-1,i-1) )
510             rkmjpim = r(k-1,j+1,i-1) + IBITS( flags(k-1,j+1,i-1), 6, 1 ) * &
511                                        ( r(k,j,i) - r(k-1,j+1,i-1) )
512             rkmjmip = r(k-1,j-1,i+1) + IBITS( flags(k-1,j-1,i+1), 6, 1 ) * &
513                                        ( r(k,j,i) - r(k-1,j-1,i+1) )
514             rkmjpip = r(k-1,j+1,i+1) + IBITS( flags(k-1,j+1,i+1), 6, 1 ) * &
515                                        ( r(k,j,i) - r(k-1,j+1,i+1) )
516
517             f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                         &
518                              8.0_wp * r(k,j,i)                            &
519                            + 4.0_wp * ( rkjim   + rkjip   +               &
520                                         rkjpi   + rkjmi   )               &
521                            + 2.0_wp * ( rkjmim  + rkjpim  +               &
522                                         rkjmip  + rkjpip  )               &
523                            + 4.0_wp * rkmji                               &
524                            + 2.0_wp * ( rkmjim  + rkmjim  +               &
525                                         rkmjpi  + rkmjmi  )               &
526                            +          ( rkmjmim + rkmjpim +               &
527                                         rkmjmip + rkmjpip )               &
528                            + 4.0_wp * r(k+1,j,i)                          &
529                            + 2.0_wp * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
530                                         r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
531                            +          ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
532                                         r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
533                                                 )
534
535!             f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                         &
536!                              8.0_wp * r(k,j,i)                            &
537!                            + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
538!                                         r(k,j+1,i)     + r(k,j-1,i)     ) &
539!                            + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
540!                                         r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
541!                            + 4.0_wp * r(k-1,j,i)                          &
542!                            + 2.0_wp * ( r(k-1,j,i-1)   + r(k-1,j,i+1)   + &
543!                                         r(k-1,j+1,i)   + r(k-1,j-1,i)   ) &
544!                            +          ( r(k-1,j-1,i-1) + r(k-1,j+1,i-1) + &
545!                                         r(k-1,j-1,i+1) + r(k-1,j+1,i+1) ) &
546!                            + 4.0_wp * r(k+1,j,i)                          &
547!                            + 2.0_wp * ( r(k+1,j,i-1)   + r(k+1,j,i+1)   + &
548!                                         r(k+1,j+1,i)   + r(k+1,j-1,i)   ) &
549!                            +          ( r(k+1,j-1,i-1) + r(k+1,j+1,i-1) + &
550!                                         r(k+1,j-1,i+1) + r(k+1,j+1,i+1) ) &
551!                                                )
552          ENDDO
553       ENDDO
554    ENDDO
555!$OMP END PARALLEL
556
557!
558!-- Horizontal boundary conditions
559    CALL exchange_horiz( f_mg, 1)
560
561    IF ( .NOT. bc_lr_cyc )  THEN
562       IF (inflow_l .OR. outflow_l)  f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
563       IF (inflow_r .OR. outflow_r)  f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
564    ENDIF
565
566    IF ( .NOT. bc_ns_cyc )  THEN
567       IF (inflow_n .OR. outflow_n)  f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
568       IF (inflow_s .OR. outflow_s)  f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
569    ENDIF
570
571!
572!-- Boundary conditions at bottom and top of the domain.
573!-- These points are not handled by the above loop. Points may be within
574!-- buildings, but that doesn't matter.
575    IF ( ibc_p_b == 1 )  THEN
576       f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
577    ELSE
578       f_mg(nzb,:,: ) = 0.0_wp
579    ENDIF
580
581    IF ( ibc_p_t == 1 )  THEN
582       f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
583    ELSE
584       f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
585    ENDIF
586
587
588END SUBROUTINE restrict
589
590
591!------------------------------------------------------------------------------!
592! Description:
593! ------------
594!> Interpolates the correction of the perturbation pressure
595!> to the next finer grid.
596!------------------------------------------------------------------------------!
597 SUBROUTINE prolong( p, temp )
598
599
600    USE control_parameters,                                                    &
601        ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,   &
602               inflow_n, inflow_r, inflow_s, outflow_l, outflow_n, outflow_r,  &
603               outflow_s
604
605    USE indices,                                                               &
606        ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
607
608    USE kinds
609
610    IMPLICIT NONE
611
612    INTEGER(iwp) ::  i  !<
613    INTEGER(iwp) ::  j  !<
614    INTEGER(iwp) ::  k  !<
615    INTEGER(iwp) ::  l  !<
616
617    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                           &
618                        nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,        &
619                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p  !<
620
621    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
622                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
623                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp  !<
624
625
626!
627!-- First, store elements of the coarser grid on the next finer grid
628    l = grid_level
629
630!$OMP PARALLEL PRIVATE (i,j,k)
631!$OMP DO
632    DO  i = nxl_mg(l-1), nxr_mg(l-1)
633       DO  j = nys_mg(l-1), nyn_mg(l-1)
634!CDIR NODEP
635          DO  k = nzb+1, nzt_mg(l-1)
636!
637!--          Points of the coarse grid are directly stored on the next finer
638!--          grid
639             temp(2*k-1,2*j,2*i) = p(k,j,i) 
640!
641!--          Points between two coarse-grid points
642             temp(2*k-1,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
643             temp(2*k-1,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
644             temp(2*k,2*j,2*i)     = 0.5_wp * ( p(k,j,i) + p(k+1,j,i) )
645!
646!--          Points in the center of the planes stretched by four points
647!--          of the coarse grid cube
648             temp(2*k-1,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) + &
649                                                   p(k,j+1,i) + p(k,j+1,i+1) )
650             temp(2*k,2*j,2*i+1)     = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) + &
651                                                   p(k+1,j,i) + p(k+1,j,i+1) )
652             temp(2*k,2*j+1,2*i)     = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) + &
653                                                   p(k+1,j,i) + p(k+1,j+1,i) )
654!
655!--          Points in the middle of coarse grid cube
656             temp(2*k,2*j+1,2*i+1) = 0.125_wp * ( p(k,j,i)     + p(k,j,i+1)   + &
657                                                  p(k,j+1,i)   + p(k,j+1,i+1) + &
658                                                  p(k+1,j,i)   + p(k+1,j,i+1) + &
659                                                  p(k+1,j+1,i) + p(k+1,j+1,i+1) )
660          ENDDO
661       ENDDO
662    ENDDO
663!$OMP END PARALLEL
664                         
665!
666!-- Horizontal boundary conditions
667    CALL exchange_horiz( temp, 1)
668
669    IF ( .NOT. bc_lr_cyc )  THEN
670       IF (inflow_l .OR. outflow_l)  temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
671       IF (inflow_r .OR. outflow_r)  temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
672    ENDIF
673
674    IF ( .NOT. bc_ns_cyc )  THEN
675       IF (inflow_n .OR. outflow_n)  temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
676       IF (inflow_s .OR. outflow_s)  temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
677    ENDIF
678
679!
680!-- Bottom and top boundary conditions
681    IF ( ibc_p_b == 1 )  THEN
682       temp(nzb,:,: ) = temp(nzb+1,:,:)
683    ELSE
684       temp(nzb,:,: ) = 0.0_wp
685    ENDIF
686
687    IF ( ibc_p_t == 1 )  THEN
688       temp(nzt_mg(l)+1,:,: ) = temp(nzt_mg(l),:,:)
689    ELSE
690       temp(nzt_mg(l)+1,:,: ) = 0.0_wp
691    ENDIF
692
693 
694 END SUBROUTINE prolong
695
696
697!------------------------------------------------------------------------------!
698! Description:
699! ------------
700!> Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
701!> 3D-Red-Black decomposition (GS-RB) is used.
702!------------------------------------------------------------------------------!
703 SUBROUTINE redblack( f_mg, p_mg )
704
705
706    USE arrays_3d,                                                             &
707        ONLY:  f1_mg, f2_mg, f3_mg
708
709    USE control_parameters,                                                    &
710        ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,   &
711               inflow_n, inflow_r, inflow_s, ngsrb, outflow_l, outflow_n,      &
712               outflow_r, outflow_s
713
714    USE cpulog,                                                                &
715        ONLY:  cpu_log, log_point_s
716
717    USE grid_variables,                                                        &
718        ONLY:  ddx2_mg, ddy2_mg
719
720    USE indices,                                                               &
721        ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3, wall_flags_4,  &
722               wall_flags_5, wall_flags_6, wall_flags_7, wall_flags_8,         &
723               wall_flags_9, wall_flags_10, nxl_mg, nxr_mg, nys_mg, nyn_mg,    &
724               nzb, nzt_mg
725
726    USE kinds
727
728    IMPLICIT NONE
729
730    INTEGER(iwp) :: color    !<
731    INTEGER(iwp) :: i        !<
732    INTEGER(iwp) :: ic       !<
733    INTEGER(iwp) :: j        !<
734    INTEGER(iwp) :: jc       !<
735    INTEGER(iwp) :: jj       !<
736    INTEGER(iwp) :: k        !<
737    INTEGER(iwp) :: l        !<
738    INTEGER(iwp) :: n        !<
739
740    LOGICAL :: unroll        !<
741
742    REAL(wp) ::  wall_left   !<
743    REAL(wp) ::  wall_north  !<
744    REAL(wp) ::  wall_right  !<
745    REAL(wp) ::  wall_south  !<
746    REAL(wp) ::  wall_total  !<
747    REAL(wp) ::  wall_top    !<
748
749    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
750                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,             &
751                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !<
752    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
753                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,             &
754                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !<
755
756    l = grid_level
757
758!
759!-- Choose flag array of this level
760    SELECT CASE ( l )
761       CASE ( 1 )
762          flags => wall_flags_1
763       CASE ( 2 )
764          flags => wall_flags_2
765       CASE ( 3 )
766          flags => wall_flags_3
767       CASE ( 4 )
768          flags => wall_flags_4
769       CASE ( 5 )
770          flags => wall_flags_5
771       CASE ( 6 )
772          flags => wall_flags_6
773       CASE ( 7 )
774          flags => wall_flags_7
775       CASE ( 8 )
776          flags => wall_flags_8
777       CASE ( 9 )
778          flags => wall_flags_9
779       CASE ( 10 )
780          flags => wall_flags_10
781    END SELECT
782
783    unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
784               MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
785
786    DO  n = 1, ngsrb
787       
788       DO  color = 1, 2
789
790          IF ( .NOT. unroll )  THEN
791
792             CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
793
794!
795!--          Without unrolling of loops, no cache optimization
796             DO  i = nxl_mg(l), nxr_mg(l), 2
797                DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2 
798                   DO  k = nzb+1, nzt_mg(l), 2
799!                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
800!                                ddx2_mg(l) * ( p_mg(k,j,i+1) + p_mg(k,j,i-1) ) &
801!                              + ddy2_mg(l) * ( p_mg(k,j+1,i) + p_mg(k,j-1,i) ) &
802!                              + f2_mg(k,l) * p_mg(k+1,j,i)                     &
803!                              + f3_mg(k,l) * p_mg(k-1,j,i) - f_mg(k,j,i)       &
804!                                                       )
805
806                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
807                             ddx2_mg(l) *                                      &
808                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
809                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
810                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
811                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
812                           + ddy2_mg(l) *                                      &
813                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
814                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
815                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
816                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
817                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
818                           + f3_mg(k,l) *                                      &
819                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
820                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
821                           - f_mg(k,j,i)               )
822                   ENDDO
823                ENDDO
824             ENDDO
825   
826             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
827                DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
828                   DO  k = nzb+1, nzt_mg(l), 2 
829                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
830                             ddx2_mg(l) *                                      &
831                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
832                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
833                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
834                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
835                           + ddy2_mg(l) *                                      &
836                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
837                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
838                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
839                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
840                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
841                           + f3_mg(k,l) *                                      &
842                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
843                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
844                           - f_mg(k,j,i)               )
845                   ENDDO
846                ENDDO
847             ENDDO
848 
849             DO  i = nxl_mg(l), nxr_mg(l), 2
850                DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
851                   DO  k = nzb+2, nzt_mg(l), 2
852                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
853                             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                           + 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) + 2 - color, nyn_mg(l), 2
874                   DO  k = nzb+2, nzt_mg(l), 2
875                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
876                             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                           + 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             CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
895
896          ELSE
897
898!
899!--          Loop unrolling along y, only one i loop for better cache use
900             CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
901             DO  ic = nxl_mg(l), nxr_mg(l), 2
902                DO  jc = nys_mg(l), nyn_mg(l), 4
903                   i  = ic
904                   jj = jc+2-color
905                   DO  k = nzb+1, nzt_mg(l), 2
906                      j = jj
907                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
908                             ddx2_mg(l) *                                      &
909                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
910                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
911                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
912                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
913                           + ddy2_mg(l) *                                      &
914                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
915                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
916                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
917                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
918                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
919                           + f3_mg(k,l) *                                      &
920                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
921                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
922                           - f_mg(k,j,i)               )
923                      j = jj+2
924                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
925                             ddx2_mg(l) *                                      &
926                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
927                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
928                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
929                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
930                           + ddy2_mg(l) *                                      &
931                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
932                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
933                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
934                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
935                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
936                           + f3_mg(k,l) *                                      &
937                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
938                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
939                           - f_mg(k,j,i)               )
940                   ENDDO
941   
942                   i  = ic+1
943                   jj = jc+color-1
944                   DO  k = nzb+1, nzt_mg(l), 2 
945                      j =jj
946                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
947                             ddx2_mg(l) *                                      &
948                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
949                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
950                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
951                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
952                           + ddy2_mg(l) *                                      &
953                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
954                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
955                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
956                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
957                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
958                           + f3_mg(k,l) *                                      &
959                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
960                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
961                           - f_mg(k,j,i)               )
962                      j = jj+2
963                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
964                             ddx2_mg(l) *                                      &
965                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
966                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
967                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
968                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
969                           + ddy2_mg(l) *                                      &
970                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
971                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
972                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
973                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
974                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
975                           + f3_mg(k,l) *                                      &
976                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
977                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
978                           - f_mg(k,j,i)               )
979                   ENDDO
980
981                   i  = ic
982                   jj = jc+color-1
983                   DO  k = nzb+2, nzt_mg(l), 2
984                      j =jj
985                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
986                             ddx2_mg(l) *                                      &
987                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
988                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
989                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
990                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
991                           + ddy2_mg(l) *                                      &
992                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
993                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
994                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
995                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
996                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
997                           + f3_mg(k,l) *                                      &
998                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
999                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1000                           - f_mg(k,j,i)               )
1001                      j = jj+2
1002                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
1003                             ddx2_mg(l) *                                      &
1004                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1005                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1006                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1007                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1008                           + ddy2_mg(l) *                                      &
1009                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1010                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1011                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1012                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1013                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
1014                           + f3_mg(k,l) *                                      &
1015                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1016                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1017                           - f_mg(k,j,i)               )
1018                   ENDDO
1019
1020                   i  = ic+1
1021                   jj = jc+2-color
1022                   DO  k = nzb+2, nzt_mg(l), 2
1023                      j =jj
1024                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
1025                             ddx2_mg(l) *                                      &
1026                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1027                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1028                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1029                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1030                           + ddy2_mg(l) *                                      &
1031                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1032                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1033                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1034                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1035                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
1036                           + f3_mg(k,l) *                                      &
1037                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1038                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1039                           - f_mg(k,j,i)               )
1040                      j = jj+2
1041                      p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (                    &
1042                             ddx2_mg(l) *                                      &
1043                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1044                                             ( p_mg(k,j,i) - p_mg(k,j,i+1) ) + &
1045                                 p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) * &
1046                                             ( p_mg(k,j,i) - p_mg(k,j,i-1) ) ) &
1047                           + ddy2_mg(l) *                                      &
1048                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1049                                             ( p_mg(k,j,i) - p_mg(k,j+1,i) ) + &
1050                                 p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) * &
1051                                             ( p_mg(k,j,i) - p_mg(k,j-1,i) ) ) &
1052                           + f2_mg(k,l) * p_mg(k+1,j,i)                        &
1053                           + f3_mg(k,l) *                                      &
1054                               ( p_mg(k-1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1055                                             ( p_mg(k,j,i) - p_mg(k-1,j,i) ) ) &
1056                           - f_mg(k,j,i)               )
1057                   ENDDO
1058
1059                ENDDO
1060             ENDDO
1061             CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
1062
1063          ENDIF
1064
1065!
1066!--       Horizontal boundary conditions
1067          CALL exchange_horiz( p_mg, 1 )
1068
1069          IF ( .NOT. bc_lr_cyc )  THEN
1070             IF ( inflow_l .OR. outflow_l )  THEN
1071                p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
1072             ENDIF
1073             IF ( inflow_r .OR. outflow_r )  THEN
1074                p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
1075             ENDIF
1076          ENDIF
1077
1078          IF ( .NOT. bc_ns_cyc )  THEN
1079             IF ( inflow_n .OR. outflow_n )  THEN
1080                p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
1081             ENDIF
1082             IF ( inflow_s .OR. outflow_s )  THEN
1083                p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
1084             ENDIF
1085          ENDIF
1086
1087!
1088!--       Bottom and top boundary conditions
1089          IF ( ibc_p_b == 1 )  THEN
1090             p_mg(nzb,:,: ) = p_mg(nzb+1,:,:)
1091          ELSE
1092             p_mg(nzb,:,: ) = 0.0_wp
1093          ENDIF
1094
1095          IF ( ibc_p_t == 1 )  THEN
1096             p_mg(nzt_mg(l)+1,:,: ) = p_mg(nzt_mg(l),:,:)
1097          ELSE
1098             p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
1099          ENDIF
1100
1101       ENDDO
1102
1103    ENDDO
1104
1105!
1106!-- Set pressure within topography and at the topography surfaces
1107!$OMP PARALLEL PRIVATE (i,j,k,wall_left,wall_north,wall_right,wall_south,wall_top,wall_total)
1108!$OMP DO
1109    DO  i = nxl_mg(l), nxr_mg(l)
1110       DO  j = nys_mg(l), nyn_mg(l) 
1111          DO  k = nzb, nzt_mg(l)
1112!
1113!--          First, set pressure inside topography to zero
1114             p_mg(k,j,i) = p_mg(k,j,i) * ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1115!
1116!--          Second, determine if the gridpoint inside topography is adjacent
1117!--          to a wall and set its value to a value given by the average of
1118!--          those values obtained from Neumann boundary condition
1119             wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
1120             wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1121             wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1122             wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1123             wall_top   = IBITS( flags(k+1,j,i), 0, 1 )
1124             wall_total = wall_left + wall_right + wall_south + wall_north + &
1125                          wall_top
1126
1127             IF ( wall_total > 0.0_wp )  THEN
1128                p_mg(k,j,i) = 1.0_wp / wall_total *                 &
1129                                     ( wall_left  * p_mg(k,j,i-1) + &
1130                                       wall_right * p_mg(k,j,i+1) + &
1131                                       wall_south * p_mg(k,j-1,i) + &
1132                                       wall_north * p_mg(k,j+1,i) + &
1133                                       wall_top   * p_mg(k+1,j,i) )
1134             ENDIF
1135          ENDDO
1136       ENDDO
1137    ENDDO
1138!$OMP END PARALLEL
1139
1140!
1141!-- One more time horizontal boundary conditions
1142    CALL exchange_horiz( p_mg, 1)
1143
1144
1145 END SUBROUTINE redblack
1146
1147
1148
1149!------------------------------------------------------------------------------!
1150! Description:
1151! ------------
1152!> Gather subdomain data from all PEs.
1153!------------------------------------------------------------------------------!
1154 SUBROUTINE mg_gather( f2, f2_sub )
1155
1156    USE control_parameters,                                                    &
1157        ONLY:  grid_level
1158
1159    USE cpulog,                                                                &
1160        ONLY:  cpu_log, log_point_s
1161
1162    USE indices,                                                               &
1163        ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1164
1165    USE kinds
1166
1167    USE pegrid
1168
1169    IMPLICIT NONE
1170
1171    INTEGER(iwp) ::  i       !<
1172    INTEGER(iwp) ::  il      !<
1173    INTEGER(iwp) ::  ir      !<
1174    INTEGER(iwp) ::  j       !<
1175    INTEGER(iwp) ::  jn      !<
1176    INTEGER(iwp) ::  js      !<
1177    INTEGER(iwp) ::  k       !<
1178    INTEGER(iwp) ::  nwords  !<
1179
1180    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
1181                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
1182                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2    !<
1183    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
1184                    nys_mg(grid_level)-1:nyn_mg(grid_level)+1,                 &
1185                    nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2_l  !<
1186
1187    REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                              &
1188                        mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,             &
1189                        mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub  !<
1190
1191
1192#if defined( __parallel )
1193    CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
1194
1195    f2_l = 0.0_wp
1196
1197!
1198!-- Store the local subdomain array on the total array
1199    js = mg_loc_ind(3,myid)
1200    IF ( south_border_pe )  js = js - 1
1201    jn = mg_loc_ind(4,myid)
1202    IF ( north_border_pe )  jn = jn + 1
1203    il = mg_loc_ind(1,myid)
1204    IF ( left_border_pe )   il = il - 1
1205    ir = mg_loc_ind(2,myid)
1206    IF ( right_border_pe )  ir = ir + 1
1207    DO  i = il, ir
1208       DO  j = js, jn
1209          DO  k = nzb, nzt_mg(grid_level)+1
1210             f2_l(k,j,i) = f2_sub(k,j,i)
1211          ENDDO
1212       ENDDO
1213    ENDDO
1214
1215!
1216!-- Find out the number of array elements of the total array
1217    nwords = SIZE( f2 )
1218
1219!
1220!-- Gather subdomain data from all PEs
1221    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1222    CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1223                        f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),   &
1224                        nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1225
1226    CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
1227#endif
1228   
1229 END SUBROUTINE mg_gather
1230
1231
1232
1233!------------------------------------------------------------------------------!
1234! Description:
1235! ------------
1236!> @todo It may be possible to improve the speed of this routine by using
1237!>       non-blocking communication
1238!------------------------------------------------------------------------------!
1239 SUBROUTINE mg_scatter( p2, p2_sub )
1240
1241    USE control_parameters,                                                    &
1242        ONLY:  grid_level
1243
1244    USE cpulog,                                                                &
1245        ONLY:  cpu_log, log_point_s
1246
1247    USE indices,                                                               &
1248        ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1249
1250    USE kinds
1251
1252    USE pegrid
1253
1254    IMPLICIT NONE
1255
1256    INTEGER(iwp) ::  nwords  !<
1257
1258    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                            &
1259                        nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,         &
1260                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
1261
1262    REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                              &
1263                        mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,             &
1264                        mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub  !<
1265
1266!
1267!-- Find out the number of array elements of the subdomain array
1268    nwords = SIZE( p2_sub )
1269
1270#if defined( __parallel )
1271    CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
1272
1273    p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1274                  mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1275
1276    CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
1277#endif
1278   
1279 END SUBROUTINE mg_scatter
1280
1281
1282!------------------------------------------------------------------------------!
1283! Description:
1284! ------------
1285!> This is where the multigrid technique takes place. V- and W- Cycle are
1286!> implemented and steered by the parameter "gamma". Parameter "nue" determines
1287!> the convergence of the multigrid iterative solution. There are nue times
1288!> RB-GS iterations. It should be set to "1" or "2", considering the time effort
1289!> one would like to invest. Last choice shows a very good converging factor,
1290!> but leads to an increase in computing time.
1291!------------------------------------------------------------------------------!
1292 RECURSIVE SUBROUTINE next_mg_level( f_mg, p_mg, p3, r )
1293
1294    USE control_parameters,                                                    &
1295        ONLY:  bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir,         &
1296               gamma_mg, grid_level, grid_level_count, ibc_p_b, ibc_p_t,       &
1297               inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level,     &
1298               mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb, outflow_l,     &
1299               outflow_n, outflow_r, outflow_s
1300
1301
1302    USE indices,                                                               &
1303        ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,         &
1304               nyn_mg, nzb, nzt, nzt_mg
1305
1306    USE kinds
1307
1308    USE pegrid
1309
1310    IMPLICIT NONE
1311
1312    INTEGER(iwp) ::  i            !<
1313    INTEGER(iwp) ::  j            !<
1314    INTEGER(iwp) ::  k            !<
1315    INTEGER(iwp) ::  nxl_mg_save  !<
1316    INTEGER(iwp) ::  nxr_mg_save  !<
1317    INTEGER(iwp) ::  nyn_mg_save  !<
1318    INTEGER(iwp) ::  nys_mg_save  !<
1319    INTEGER(iwp) ::  nzt_mg_save  !<
1320
1321    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
1322                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,             &
1323                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg  !<
1324    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
1325                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,             &
1326                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg  !<
1327    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
1328                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,             &
1329                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3    !<
1330    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                              &
1331                        nys_mg(grid_level)-1:nyn_mg(grid_level)+1,             &
1332                        nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r     !<
1333
1334    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                            &
1335                        nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,         &
1336                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2  !<
1337    REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                            &
1338                        nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,         &
1339                        nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !<
1340
1341    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub  !<
1342    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p2_sub  !<
1343
1344!
1345!-- Restriction to the coarsest grid
1346 10 IF ( grid_level == 1 )  THEN
1347
1348!
1349!--    Solution on the coarsest grid. Double the number of Gauss-Seidel
1350!--    iterations in order to get a more accurate solution.
1351       ngsrb = 2 * ngsrb
1352
1353       CALL redblack( f_mg, p_mg )
1354
1355       ngsrb = ngsrb / 2
1356
1357
1358    ELSEIF ( grid_level /= 1 )  THEN
1359
1360       grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1361
1362!
1363!--    Solution on the actual grid level
1364       CALL redblack( f_mg, p_mg )
1365
1366!
1367!--    Determination of the actual residual
1368       CALL resid( f_mg, p_mg, r )
1369
1370!
1371!--    Restriction of the residual (finer grid values!) to the next coarser
1372!--    grid. Therefore, the grid level has to be decremented now. nxl..nzt have
1373!--    to be set to the coarse grid values, because these variables are needed
1374!--    for the exchange of ghost points in routine exchange_horiz
1375       grid_level = grid_level - 1
1376       nxl = nxl_mg(grid_level)
1377       nys = nys_mg(grid_level)
1378       nxr = nxr_mg(grid_level)
1379       nyn = nyn_mg(grid_level)
1380       nzt = nzt_mg(grid_level)
1381
1382       IF ( grid_level == mg_switch_to_pe0_level )  THEN
1383
1384!
1385!--       From this level on, calculations are done on PE0 only.
1386!--       First, carry out restriction on the subdomain.
1387!--       Therefore, indices of the level have to be changed to subdomain values
1388!--       in between (otherwise, the restrict routine would expect
1389!--       the gathered array)
1390
1391          nxl_mg_save = nxl_mg(grid_level)
1392          nxr_mg_save = nxr_mg(grid_level)
1393          nys_mg_save = nys_mg(grid_level)
1394          nyn_mg_save = nyn_mg(grid_level)
1395          nzt_mg_save = nzt_mg(grid_level)
1396          nxl_mg(grid_level) = mg_loc_ind(1,myid)
1397          nxr_mg(grid_level) = mg_loc_ind(2,myid)
1398          nys_mg(grid_level) = mg_loc_ind(3,myid)
1399          nyn_mg(grid_level) = mg_loc_ind(4,myid)
1400          nzt_mg(grid_level) = mg_loc_ind(5,myid)
1401          nxl = mg_loc_ind(1,myid)
1402          nxr = mg_loc_ind(2,myid)
1403          nys = mg_loc_ind(3,myid)
1404          nyn = mg_loc_ind(4,myid)
1405          nzt = mg_loc_ind(5,myid)
1406
1407          ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
1408                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
1409                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
1410
1411          CALL restrict( f2_sub, r )
1412
1413!
1414!--       Restore the correct indices of this level
1415          nxl_mg(grid_level) = nxl_mg_save
1416          nxr_mg(grid_level) = nxr_mg_save
1417          nys_mg(grid_level) = nys_mg_save
1418          nyn_mg(grid_level) = nyn_mg_save
1419          nzt_mg(grid_level) = nzt_mg_save
1420          nxl = nxl_mg(grid_level)
1421          nxr = nxr_mg(grid_level)
1422          nys = nys_mg(grid_level)
1423          nyn = nyn_mg(grid_level)
1424          nzt = nzt_mg(grid_level)
1425!
1426!--       Gather all arrays from the subdomains on PE0
1427          CALL mg_gather( f2, f2_sub )
1428
1429!
1430!--       Set switch for routine exchange_horiz, that no ghostpoint exchange
1431!--       has to be carried out from now on
1432          mg_switch_to_pe0 = .TRUE.
1433
1434!
1435!--       In case of non-cyclic lateral boundary conditions, both in- and
1436!--       outflow conditions have to be used on all PEs after the switch,
1437!--       because then they have the total domain.
1438          IF ( bc_lr_dirrad )  THEN
1439             inflow_l  = .TRUE.
1440             inflow_r  = .FALSE.
1441             outflow_l = .FALSE.
1442             outflow_r = .TRUE.
1443          ELSEIF ( bc_lr_raddir )  THEN
1444             inflow_l  = .FALSE.
1445             inflow_r  = .TRUE.
1446             outflow_l = .TRUE.
1447             outflow_r = .FALSE.
1448          ENDIF
1449
1450          IF ( bc_ns_dirrad )  THEN
1451             inflow_n  = .TRUE.
1452             inflow_s  = .FALSE.
1453             outflow_n = .FALSE.
1454             outflow_s = .TRUE.
1455          ELSEIF ( bc_ns_raddir )  THEN
1456             inflow_n  = .FALSE.
1457             inflow_s  = .TRUE.
1458             outflow_n = .TRUE.
1459             outflow_s = .FALSE.
1460          ENDIF
1461
1462          DEALLOCATE( f2_sub )
1463
1464       ELSE
1465
1466          CALL restrict( f2, r )
1467
1468       ENDIF
1469
1470       p2 = 0.0_wp
1471
1472!
1473!--    Repeat the same procedure till the coarsest grid is reached
1474       CALL next_mg_level( f2, p2, p3, r )
1475
1476    ENDIF
1477
1478!
1479!-- Now follows the prolongation
1480    IF ( grid_level >= 2 )  THEN
1481
1482!
1483!--    Prolongation of the new residual. The values are transferred
1484!--    from the coarse to the next finer grid.
1485       IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
1486
1487#if defined( __parallel )
1488!
1489!--       At this level, the new residual first has to be scattered from
1490!--       PE0 to the other PEs
1491          ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
1492                    mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
1493                    mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
1494
1495          CALL mg_scatter( p2, p2_sub )
1496
1497!
1498!--       Therefore, indices of the previous level have to be changed to
1499!--       subdomain values in between (otherwise, the prolong routine would
1500!--       expect the gathered array)
1501          nxl_mg_save = nxl_mg(grid_level-1)
1502          nxr_mg_save = nxr_mg(grid_level-1)
1503          nys_mg_save = nys_mg(grid_level-1)
1504          nyn_mg_save = nyn_mg(grid_level-1)
1505          nzt_mg_save = nzt_mg(grid_level-1)
1506          nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
1507          nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
1508          nys_mg(grid_level-1) = mg_loc_ind(3,myid)
1509          nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
1510          nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
1511
1512!
1513!--       Set switch for routine exchange_horiz, that ghostpoint exchange
1514!--       has to be carried again out from now on
1515          mg_switch_to_pe0 = .FALSE.
1516
1517!
1518!--       For non-cyclic lateral boundary conditions, restore the
1519!--       in-/outflow conditions
1520          inflow_l  = .FALSE.;  inflow_r  = .FALSE.
1521          inflow_n  = .FALSE.;  inflow_s  = .FALSE.
1522          outflow_l = .FALSE.;  outflow_r = .FALSE.
1523          outflow_n = .FALSE.;  outflow_s = .FALSE.
1524
1525          IF ( pleft == MPI_PROC_NULL )  THEN
1526             IF ( bc_lr_dirrad )  THEN
1527                inflow_l  = .TRUE.
1528             ELSEIF ( bc_lr_raddir )  THEN
1529                outflow_l = .TRUE.
1530             ENDIF
1531          ENDIF
1532
1533          IF ( pright == MPI_PROC_NULL )  THEN
1534             IF ( bc_lr_dirrad )  THEN
1535                outflow_r = .TRUE.
1536             ELSEIF ( bc_lr_raddir )  THEN
1537                inflow_r  = .TRUE.
1538             ENDIF
1539          ENDIF
1540
1541          IF ( psouth == MPI_PROC_NULL )  THEN
1542             IF ( bc_ns_dirrad )  THEN
1543                outflow_s = .TRUE.
1544             ELSEIF ( bc_ns_raddir )  THEN
1545                inflow_s  = .TRUE.
1546             ENDIF
1547          ENDIF
1548
1549          IF ( pnorth == MPI_PROC_NULL )  THEN
1550             IF ( bc_ns_dirrad )  THEN
1551                inflow_n  = .TRUE.
1552             ELSEIF ( bc_ns_raddir )  THEN
1553                outflow_n = .TRUE.
1554             ENDIF
1555          ENDIF
1556
1557          CALL prolong( p2_sub, p3 )
1558
1559!
1560!--       Restore the correct indices of the previous level
1561          nxl_mg(grid_level-1) = nxl_mg_save
1562          nxr_mg(grid_level-1) = nxr_mg_save
1563          nys_mg(grid_level-1) = nys_mg_save
1564          nyn_mg(grid_level-1) = nyn_mg_save
1565          nzt_mg(grid_level-1) = nzt_mg_save
1566
1567          DEALLOCATE( p2_sub )
1568#endif
1569
1570       ELSE
1571
1572          CALL prolong( p2, p3 )
1573
1574       ENDIF
1575
1576!
1577!--    Computation of the new pressure correction. Therefore,
1578!--    values from prior grids are added up automatically stage by stage.
1579       DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
1580          DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
1581             DO  k = nzb, nzt_mg(grid_level)+1 
1582                p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
1583             ENDDO
1584          ENDDO
1585       ENDDO
1586
1587!
1588!--    Relaxation of the new solution
1589       CALL redblack( f_mg, p_mg )
1590
1591    ENDIF
1592
1593
1594!
1595!-- The following few lines serve the steering of the multigrid scheme
1596    IF ( grid_level == maximum_grid_level )  THEN
1597
1598       GOTO 20
1599
1600    ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
1601             grid_level_count(grid_level) /= gamma_mg )  THEN
1602
1603       GOTO 10
1604
1605    ENDIF
1606
1607!
1608!-- Reset counter for the next call of poismg
1609    grid_level_count(grid_level) = 0
1610
1611!
1612!-- Continue with the next finer level. nxl..nzt have to be
1613!-- set to the finer grid values, because these variables are needed for the
1614!-- exchange of ghost points in routine exchange_horiz
1615    grid_level = grid_level + 1
1616    nxl = nxl_mg(grid_level)
1617    nxr = nxr_mg(grid_level)
1618    nys = nys_mg(grid_level)
1619    nyn = nyn_mg(grid_level)
1620    nzt = nzt_mg(grid_level)
1621
1622 20 CONTINUE
1623
1624 END SUBROUTINE next_mg_level
Note: See TracBrowser for help on using the repository browser.