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

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

last commit documented

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