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

Last change on this file since 2339 was 2233, checked in by suehring, 7 years ago

last commit documented

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