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

Last change on this file since 2000 was 2000, checked in by knoop, 5 years ago

Forced header and separation lines into 80 columns

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