source: palm/tags/release-4.0/SOURCE/poismg.f90 @ 3800

Last change on this file since 3800 was 1354, checked in by heinze, 10 years ago

last commit documented

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