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

Last change on this file since 1834 was 1818, checked in by maronga, 5 years ago

last commit documented / copyright update

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