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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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