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

Last change on this file since 1322 was 1322, checked in by raasch, 8 years ago

REAL functions and a lot of REAL constants provided with KIND-attribute,
some small bugfixes

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