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

Last change on this file since 2232 was 2232, checked in by suehring, 4 years ago

Adjustments according new topography and surface-modelling concept implemented

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