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

Last change on this file since 2012 was 2001, checked in by knoop, 8 years ago

last commit documented

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