source: palm/trunk/SOURCE/poismg_noopt.f90 @ 1985

Last change on this file since 1985 was 1935, checked in by hellstea, 8 years ago

last commit documented

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