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

Last change on this file since 1931 was 1931, checked in by suehring, 8 years ago

Rename multigrid into multigrid_noopt and multigrid_fast into multigrid, subroutines poismg is renamed into poismg_noopt and poismg_fast_mod into poismg_mod

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