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

Last change on this file since 2101 was 2101, checked in by suehring, 5 years ago

last commit documented

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