source: palm/trunk/SOURCE/poismg_fast.f90 @ 1611

Last change on this file since 1611 was 1610, checked in by maronga, 9 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 113.2 KB
Line 
1 MODULE poismg_mod
2
3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later 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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Attention: Loop unrolling and cache optimization in SOR-Red/Black method
21!            still does not give the expected speedup! Further work required.
22!
23! Current revisions:
24! -----------------
25!
26!
27! Former revisions:
28! -----------------
29! $Id: poismg_fast.f90 1610 2015-07-03 15:38:36Z maronga $
30!
31! 1609 2015-07-03 15:37:58Z maronga
32! Bugfix: allow compilation without __parallel.
33!
34! 1575 2015-03-27 09:56:27Z raasch
35! Initial revision.
36! Routine re-written and optimised based on poismg.
37!
38! Following optimisations have been made:
39! - vectorisation (for Intel-CPUs) of the red-black algorithm by resorting
40!   array elements with even and odd indices
41! - explicit boundary conditions for building walls removed (solver is
42!   running through the buildings
43! - reduced data transfer in case of ghost point exchange, because only
44!   "red" or "black" data points need to be exchanged. This is not applied
45!   for coarser grid levels, since for then the transfer time is latency bound
46!
47!
48! Description:
49! ------------
50! Solves the Poisson equation for the perturbation pressure with a multigrid
51! V- or W-Cycle scheme.
52!
53! This multigrid method was originally developed for PALM by Joerg Uhlenbrock,
54! September 2000 - July 2001. It has been optimised for speed by Klaus
55! Ketelsen in November 2014.
56!------------------------------------------------------------------------------!
57
58    USE cpulog,                                                                &
59        ONLY:  cpu_log, log_point_s
60
61    USE kinds
62
63    USE pegrid
64
65    PRIVATE
66
67    INTEGER, SAVE                             ::  ind_even_odd    !: border index between even and odd k index
68    INTEGER, DIMENSION(:), SAVE, ALLOCATABLE  ::  even_odd_level  !: stores ind_even_odd for all MG levels
69
70    REAL(wp), DIMENSION(:,:), SAVE, ALLOCATABLE ::  f1_mg_b, f2_mg_b, f3_mg_b  !: blocked version of f1_mg ...
71
72    INTERFACE poismg_fast
73       MODULE PROCEDURE poismg_fast
74    END INTERFACE poismg_fast
75
76    INTERFACE sort_k_to_even_odd_blocks
77       MODULE PROCEDURE sort_k_to_even_odd_blocks
78       MODULE PROCEDURE sort_k_to_even_odd_blocks_int
79       MODULE PROCEDURE sort_k_to_even_odd_blocks_1d
80    END INTERFACE sort_k_to_even_odd_blocks
81
82    PUBLIC poismg_fast
83
84 CONTAINS
85
86    SUBROUTINE poismg_fast( r )
87
88       USE arrays_3d,                                                          &
89           ONLY:  d, p_loc
90
91       USE control_parameters,                                                 &
92           ONLY:  gathered_size, grid_level, grid_level_count,                 &
93                  maximum_grid_level, message_string, mgcycles, mg_cycles,     &
94                  mg_switch_to_pe0_level, residual_limit, subdomain_size
95
96       USE cpulog,                                                             &
97           ONLY:  cpu_log, log_point_s
98
99       USE indices,                                                            &
100           ONLY:  nxl, nxlg, nxl_mg, nxr, nxrg, nxr_mg, nys, nysg, nys_mg, nyn,&
101                  nyng, nyn_mg, nzb, nzt, nzt_mg
102
103       IMPLICIT NONE
104
105       REAL(wp) ::  maxerror          !:
106       REAL(wp) ::  maximum_mgcycles  !:
107       REAL(wp) ::  residual_norm     !:
108
109       REAL(wp), DIMENSION(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ::  r  !:
110
111       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p3  !:
112
113
114       CALL cpu_log( log_point_s(29), 'poismg_fast', 'start' )
115!
116!--    Initialize arrays and variables used in this subroutine
117
118!--    If the number of grid points of the gathered grid, which is collected
119!--    on PE0, is larger than the number of grid points of an PE, than array
120!--    p3 will be enlarged.
121       IF ( gathered_size > subdomain_size )  THEN
122          ALLOCATE( p3(nzb:nzt_mg(mg_switch_to_pe0_level)+1,nys_mg(            &
123                    mg_switch_to_pe0_level)-1:nyn_mg(mg_switch_to_pe0_level)+1,&
124                    nxl_mg(mg_switch_to_pe0_level)-1:nxr_mg(                   &
125                    mg_switch_to_pe0_level)+1) )
126       ELSE
127          ALLOCATE ( p3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
128       ENDIF
129
130       p3 = 0.0_wp
131
132 
133!
134!--    Ghost boundaries have to be added to divergence array.
135!--    Exchange routine needs to know the grid level!
136       grid_level = maximum_grid_level
137       CALL exchange_horiz( d, 1)
138       d(nzb,:,:) = d(nzb+1,:,:)
139
140!
141!--    Initiation of the multigrid scheme. Does n cycles until the
142!--    residual is smaller than the given limit. The accuracy of the solution
143!--    of the poisson equation will increase with the number of cycles.
144!--    If the number of cycles is preset by the user, this number will be
145!--    carried out regardless of the accuracy.
146       grid_level_count =  0
147       mgcycles         =  0
148       IF ( mg_cycles == -1 )  THEN
149          maximum_mgcycles = 0
150          residual_norm    = 1.0_wp
151       ELSE
152          maximum_mgcycles = mg_cycles
153          residual_norm    = 0.0_wp
154       ENDIF
155
156!
157!--    Initial settings for sorting k-dimension from sequential order (alternate
158!--    even/odd) into blocks of even and odd or vice versa
159       CALL init_even_odd_blocks
160
161!
162!--    Sort input arrays in even/odd blocks along k-dimension
163       CALL sort_k_to_even_odd_blocks( d, grid_level )
164       CALL sort_k_to_even_odd_blocks( p_loc, grid_level )
165
166!
167!--    The complete multigrid cycles are running in block mode, i.e. over
168!--    seperate data blocks of even and odd indices
169       DO WHILE ( residual_norm > residual_limit  .OR. &
170                  mgcycles < maximum_mgcycles )
171 
172          CALL next_mg_level_fast( d, p_loc, p3, r)
173
174!
175!--       Calculate the residual if the user has not preset the number of
176!--       cycles to be performed
177          IF ( maximum_mgcycles == 0 )  THEN
178             CALL resid_fast( d, p_loc, r )
179             maxerror = SUM( r(nzb+1:nzt,nys:nyn,nxl:nxr)**2 )
180
181#if defined( __parallel )
182             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
183                CALL MPI_ALLREDUCE( maxerror, residual_norm, 1, MPI_REAL,      &
184                                    MPI_SUM, comm2d, ierr)
185#else
186                residual_norm = maxerror
187#endif
188             residual_norm = SQRT( residual_norm )
189          ENDIF
190
191          mgcycles = mgcycles + 1
192
193!
194!--       If the user has not limited the number of cycles, stop the run in case
195!--       of insufficient convergence
196          IF ( mgcycles > 1000  .AND.  mg_cycles == -1 )  THEN
197             message_string = 'no sufficient convergence within 1000 cycles'
198             CALL message( 'poismg_fast', 'PA0283', 1, 2, 0, 6, 0 )
199          ENDIF
200
201       ENDDO
202
203       DEALLOCATE( p3 )
204!
205!--    Result has to be sorted back from even/odd blocks to sequential order
206       CALL sort_k_to_sequential( p_loc )
207!
208!--    Unset the grid level. Variable is used to determine the MPI datatypes for
209!--    ghost point exchange
210       grid_level = 0
211
212       CALL cpu_log( log_point_s(29), 'poismg_fast', 'stop' )
213
214    END SUBROUTINE poismg_fast
215
216
217
218    SUBROUTINE resid_fast( f_mg, p_mg, r )
219
220!------------------------------------------------------------------------------!
221! Description:
222! ------------
223! Computes the residual of the perturbation pressure.
224!------------------------------------------------------------------------------!
225
226       USE arrays_3d,                                                          &
227           ONLY:  f1_mg, f2_mg, f3_mg
228
229       USE control_parameters,                                                 &
230           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
231                  inflow_n, inflow_r, inflow_s, masking_method, outflow_l,     &
232                  outflow_n, outflow_r, outflow_s, topography
233
234       USE grid_variables,                                                     &
235           ONLY:  ddx2_mg, ddy2_mg
236
237       USE indices,                                                            &
238           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
239                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
240                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
241                  nys_mg, nyn_mg, nzb, nzt_mg
242
243       IMPLICIT NONE
244
245       INTEGER(iwp) ::  i
246       INTEGER(iwp) ::  j
247       INTEGER(iwp) ::  k
248       INTEGER(iwp) ::  l
249       INTEGER(iwp) :: km1      !:
250       INTEGER(iwp) :: kp1      !:
251
252       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
253                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
254                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !:
255       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
256                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
257                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
258       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
259                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
260                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  r     !:
261
262!
263!--    Calculate the residual
264       l = grid_level
265
266       CALL cpu_log( log_point_s(53), 'resid', 'start' )
267       IF ( topography == 'flat'  .OR.  masking_method )  THEN
268          !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
269          !$OMP DO
270          DO  i = nxl_mg(l), nxr_mg(l)
271             DO  j = nys_mg(l), nyn_mg(l)
272                !DIR$ IVDEP
273                DO k = ind_even_odd+1, nzt_mg(l)
274                   km1 = k-ind_even_odd-1
275                   kp1 = k-ind_even_odd
276                   r(k,j,i) = f_mg(k,j,i)                                 &
277                      - ddx2_mg(l) *                                      &
278                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                 &
279                      - ddy2_mg(l) *                                      &
280                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                  &
281                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                      &
282                      - f3_mg_b(k,l) * p_mg(km1,j,i)                      &
283                      + f1_mg_b(k,l) * p_mg(k,j,i)
284                ENDDO
285                !DIR$ IVDEP
286                DO k = nzb+1, ind_even_odd
287                   km1 = k+ind_even_odd
288                   kp1 = k+ind_even_odd+1
289                   r(k,j,i) = f_mg(k,j,i)                                 &
290                      - ddx2_mg(l) *                                      &
291                      ( p_mg(k,j,i+1) +  p_mg(k,j,i-1)  )                 &
292                      - ddy2_mg(l) *                                      &
293                      ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )                  &
294                      - f2_mg_b(k,l) * p_mg(kp1,j,i)                      &
295                      - f3_mg_b(k,l) * p_mg(km1,j,i)                      &
296                      + f1_mg_b(k,l) * p_mg(k,j,i)
297                ENDDO
298             ENDDO
299          ENDDO
300          !$OMP END PARALLEL
301       ELSE
302!
303!--       Choose flag array of this level
304          SELECT CASE ( l )
305             CASE ( 1 )
306                flags => wall_flags_1
307             CASE ( 2 )
308                flags => wall_flags_2
309             CASE ( 3 )
310                flags => wall_flags_3
311             CASE ( 4 )
312                flags => wall_flags_4
313             CASE ( 5 )
314                flags => wall_flags_5
315             CASE ( 6 )
316                flags => wall_flags_6
317             CASE ( 7 )
318                flags => wall_flags_7
319             CASE ( 8 )
320                flags => wall_flags_8
321             CASE ( 9 )
322                flags => wall_flags_9
323             CASE ( 10 )
324                flags => wall_flags_10
325          END SELECT
326
327          !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
328          !$OMP DO
329          DO  i = nxl_mg(l), nxr_mg(l)
330             DO  j = nys_mg(l), nyn_mg(l)
331
332                !DIR$ IVDEP
333                DO  k = ind_even_odd+1, nzt_mg(l)
334                   km1 = k-ind_even_odd-1
335                   kp1 = k-ind_even_odd
336
337                   r(k,j,i) = f_mg(k,j,i)                                         &
338                    - ddx2_mg(l) *                                                &
339                    ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 ))   &
340                    + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
341                    - ddy2_mg(l)                                                  &
342                    * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
343                    + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
344                    - f2_mg(k,l) * p_mg(kp1,j,i)                                  &
345                    - f3_mg(k,l) *                                                &
346                      MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 ))   &
347                    + f1_mg(k,l) * p_mg(k,j,i)
348                ENDDO
349
350                !DIR$ IVDEP
351                DO  k = nzb+1, ind_even_odd
352                   km1 = k+ind_even_odd
353                   kp1 = k+ind_even_odd+1
354
355                   r(k,j,i) = f_mg(k,j,i)                                         &
356                    - ddx2_mg(l) *                                                &
357                    ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5 ))   &
358                    + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4 )) ) &
359                    - ddy2_mg(l)                                                  &
360                    * ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3 )) &
361                    + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2 )) ) &
362                    - f2_mg(k,l) * p_mg(kp1,j,i)                                  &
363                    - f3_mg(k,l) *                                                &
364                      MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0 ))   &
365                    + f1_mg(k,l) * p_mg(k,j,i)
366                ENDDO
367!
368!--             The residual within topography should be zero
369                WHERE( BTEST(flags(nzb+1:nzt_mg(l),j,i), 6) )
370                   r(nzb+1:nzt_mg(l),j,i) = 0.0
371                END WHERE
372
373             ENDDO
374          ENDDO
375          !$OMP END PARALLEL
376
377       ENDIF
378
379!
380!--    Horizontal boundary conditions
381       CALL exchange_horiz( r, 1)
382
383       IF ( .NOT. bc_lr_cyc )  THEN
384          IF ( inflow_l .OR. outflow_l )  r(:,:,nxl_mg(l)-1) = r(:,:,nxl_mg(l))
385          IF ( inflow_r .OR. outflow_r )  r(:,:,nxr_mg(l)+1) = r(:,:,nxr_mg(l))
386       ENDIF
387
388       IF ( .NOT. bc_ns_cyc )  THEN
389          IF ( inflow_n .OR. outflow_n )  r(:,nyn_mg(l)+1,:) = r(:,nyn_mg(l),:)
390          IF ( inflow_s .OR. outflow_s )  r(:,nys_mg(l)-1,:) = r(:,nys_mg(l),:)
391       ENDIF
392
393!
394!--    Boundary conditions at bottom and top of the domain.outflow_l, outflow_n,
395!--    These points are not handled by the above loop. Points may be within
396!--    buildings, but that doesn't matter.
397       IF ( ibc_p_b == 1 )  THEN
398          r(nzb,:,: ) = r(nzb+1,:,:)
399       ELSE
400          r(nzb,:,: ) = 0.0_wp
401       ENDIF
402
403       IF ( ibc_p_t == 1 )  THEN
404          r(nzt_mg(l)+1,:,: ) = r(nzt_mg(l),:,:)
405       ELSE
406          r(nzt_mg(l)+1,:,: ) = 0.0_wp
407       ENDIF
408
409       CALL cpu_log( log_point_s(53), 'resid', 'stop' )
410
411    END SUBROUTINE resid_fast
412
413
414
415    SUBROUTINE restrict_fast( f_mg, r )
416
417!------------------------------------------------------------------------------!
418! Description:
419! ------------
420! Interpolates the residual on the next coarser grid with "full weighting"
421! scheme
422!------------------------------------------------------------------------------!
423
424       USE control_parameters,                                                 &
425           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
426                  inflow_n, inflow_r, inflow_s, masking_method, outflow_l,     &
427                  outflow_n, outflow_r, outflow_s, topography
428
429       USE indices,                                                            &
430           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
431                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
432                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
433                  nys_mg, nyn_mg, nzb, nzt_mg
434
435       IMPLICIT NONE
436
437       INTEGER(iwp) ::  i    !:
438       INTEGER(iwp) ::  ic   !:
439       INTEGER(iwp) ::  j    !:
440       INTEGER(iwp) ::  jc   !:
441       INTEGER(iwp) ::  k    !:
442       INTEGER(iwp) ::  kc   !:
443       INTEGER(iwp) ::  l    !:
444       INTEGER(iwp) ::  km1  !:
445       INTEGER(iwp) ::  kp1  !:
446
447       REAL(wp) ::  rkjim    !:
448       REAL(wp) ::  rkjip    !:
449       REAL(wp) ::  rkjmi    !:
450       REAL(wp) ::  rkjmim   !:
451       REAL(wp) ::  rkjmip   !:
452       REAL(wp) ::  rkjpi    !:
453       REAL(wp) ::  rkjpim   !:
454       REAL(wp) ::  rkjpip   !:
455       REAL(wp) ::  rkmji    !:
456       REAL(wp) ::  rkmjim   !:
457       REAL(wp) ::  rkmjip   !:
458       REAL(wp) ::  rkmjmi   !:
459       REAL(wp) ::  rkmjmim  !:
460       REAL(wp) ::  rkmjmip  !:
461       REAL(wp) ::  rkmjpi   !:
462       REAL(wp) ::  rkmjpim  !:
463       REAL(wp) ::  rkmjpip  !:
464
465       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
466                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
467                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !:
468
469       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level+1)+1,                         &
470                           nys_mg(grid_level+1)-1:nyn_mg(grid_level+1)+1,      &
471                           nxl_mg(grid_level+1)-1:nxr_mg(grid_level+1)+1) ::  r !:
472
473!
474!--    Interpolate the residual
475       l = grid_level
476
477       CALL cpu_log( log_point_s(54), 'restrict', 'start' )
478
479       IF ( topography == 'flat'  .OR.  masking_method )  THEN
480!
481!--       No wall treatment
482          !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1)
483          DO  ic = nxl_mg(l), nxr_mg(l)
484             i = 2*ic
485             !$OMP DO SCHEDULE( STATIC )
486             DO  jc = nys_mg(l), nyn_mg(l)
487!
488!--             Calculation for the first point along k
489                j  = 2*jc
490                k   = ind_even_odd+1
491                kp1 = k-ind_even_odd
492                kc  = k-ind_even_odd
493                f_mg(kc,jc,ic) =  1.0_wp / 64.0_wp * (                         &
494                                  8.0_wp * r(k,j,i)                            &
495                               +  4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
496                                             r(k,j+1,i)     + r(k,j-1,i)     ) &
497                               +  2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
498                                             r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
499                               + 16.0_wp * r(k,j,i)                            &
500                               +  4.0_wp * r(kp1,j,i)                          &
501                               +  2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
502                                             r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
503                               +           ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
504                                             r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
505                                                      )
506!
507!--             Calculation for the other points along k
508                !DIR$ IVDEP
509                DO k = ind_even_odd+2, nzt_mg(l+1)    ! Fine grid at this point
510                   km1 = k-ind_even_odd-1
511                   kp1 = k-ind_even_odd
512                   kc  = k-ind_even_odd               ! Coarse grid index
513
514                   f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
515                                  8.0_wp * r(k,j,i)                            &
516                                + 4.0_wp * ( r(k,j,i-1)     + r(k,j,i+1)     + &
517                                             r(k,j+1,i)     + r(k,j-1,i)     ) &
518                                + 2.0_wp * ( r(k,j-1,i-1)   + r(k,j+1,i-1)   + &
519                                             r(k,j-1,i+1)   + r(k,j+1,i+1)   ) &
520                                + 4.0_wp * r(km1,j,i)                          &
521                                + 2.0_wp * ( r(km1,j,i-1)   + r(km1,j,i+1)   + &
522                                             r(km1,j+1,i)   + r(km1,j-1,i)   ) &
523                                +          ( r(km1,j-1,i-1) + r(km1,j+1,i-1) + &
524                                             r(km1,j-1,i+1) + r(km1,j+1,i+1) ) &
525                                + 4.0_wp * r(kp1,j,i)                          &
526                                + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   + &
527                                             r(kp1,j+1,i)   + r(kp1,j-1,i)   ) &
528                                +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) + &
529                                             r(kp1,j-1,i+1) + r(kp1,j+1,i+1) ) &
530                                           )
531                ENDDO
532             ENDDO
533             !$OMP ENDDO nowait
534          ENDDO
535          !$OMP END PARALLEL
536
537       ELSE
538!
539!--       Choose flag array of the upper level
540          SELECT CASE ( l+1 )
541             CASE ( 1 )
542                flags => wall_flags_1
543             CASE ( 2 )
544                flags => wall_flags_2
545             CASE ( 3 )
546                flags => wall_flags_3
547             CASE ( 4 )
548                flags => wall_flags_4
549             CASE ( 5 )
550                flags => wall_flags_5
551             CASE ( 6 )
552                flags => wall_flags_6
553             CASE ( 7 )
554                flags => wall_flags_7
555             CASE ( 8 )
556                flags => wall_flags_8
557             CASE ( 9 )
558                flags => wall_flags_9
559             CASE ( 10 )
560                flags => wall_flags_10
561          END SELECT
562   
563          !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,kc,km1,kp1,rkjim,rkjip,rkjpi,    &
564          !$OMP                   rkjmim,rkjpim,rkjmip,rkjpip,rkmji,rkmjim,    &
565          !$OMP                   rkmjip,rkmjpi,rkmjmi,rkmjmim,rkmjpim,        &
566          !$OMP                   rkmjmip,rkmjpip)
567          !$OMP DO
568          DO  ic = nxl_mg(l), nxr_mg(l)
569             i = 2*ic
570             DO  jc = nys_mg(l), nyn_mg(l)
571                j = 2*jc
572                !DIR$ IVDEP
573                DO k = ind_even_odd+1, nzt_mg(l+1)      !fine grid at this point
574
575                    km1 = k-ind_even_odd-1
576                    kp1 = k-ind_even_odd
577                    kc  = k-ind_even_odd                !Coarse grid index
578!
579!--                 Use implicit Neumann BCs if the respective gridpoint is
580!--                 inside the building
581                    rkjim   = MERGE(r(k,j,i),r(k,j,i-1),BTEST( flags(k,j,i-1), 6))
582                    rkjip   = MERGE(r(k,j,i),r(k,j,i+1),BTEST( flags(k,j,i+1), 6))
583                    rkjpi   = MERGE(r(k,j,i),r(k,j+1,i),BTEST( flags(k,j+1,i), 6))
584                    rkjmi   = MERGE(r(k,j,i),r(k,j-1,i),BTEST( flags(k,j-1,i), 6))
585
586                    rkjmim  = MERGE(r(k,j,i),r(k,j-1,i-1),BTEST( flags(k,j-1,i-1), 6))
587                    rkjpim  = MERGE(r(k,j,i),r(k,j+1,i-1),BTEST( flags(k,j+1,i-1), 6))
588                    rkjmip  = MERGE(r(k,j,i),r(k,j-1,i+1),BTEST( flags(k,j-1,i+1), 6))
589                    rkjpip  = MERGE(r(k,j,i),r(k,j+1,i+1),BTEST( flags(k,j+1,i+1), 6))
590
591                    rkmji   = MERGE(r(k,j,i),r(km1,j,i)  ,BTEST( flags(km1,j,i)  , 6))
592                    rkmjim  = MERGE(r(k,j,i),r(km1,j,i-1),BTEST( flags(km1,j,i-1), 6))
593                    rkmjip  = MERGE(r(k,j,i),r(km1,j,i+1),BTEST( flags(km1,j,i+1), 6))
594                    rkmjpi  = MERGE(r(k,j,i),r(km1,j+1,i),BTEST( flags(km1,j+1,i), 6))
595                    rkmjmi  = MERGE(r(k,j,i),r(km1,j-1,i),BTEST( flags(km1,j-1,i), 6))
596
597                    rkmjmim = MERGE(r(k,j,i),r(km1,j-1,i-1),BTEST( flags(km1,j-1,i-1), 6))
598                    rkmjpim = MERGE(r(k,j,i),r(km1,j+1,i-1),BTEST( flags(km1,j+1,i-1), 6))
599                    rkmjmip = MERGE(r(k,j,i),r(km1,j-1,i+1),BTEST( flags(km1,j-1,i+1), 6))
600                    rkmjpip = MERGE(r(k,j,i),r(km1,j+1,i+1),BTEST( flags(km1,j+1,i+1), 6))
601
602                    f_mg(kc,jc,ic) = 1.0_wp / 64.0_wp * (                      &
603                           8.0_wp * r(k,j,i)                                   &
604                         + 4.0_wp * ( rkjim   + rkjip   +  rkjpi   + rkjmi   ) &
605                         + 2.0_wp * ( rkjmim  + rkjpim  +  rkjmip  + rkjpip  ) &
606                         + 4.0_wp * rkmji                                      &
607                         + 2.0_wp * ( rkmjim  + rkmjip  +  rkmjpi  + rkmjmi  ) &
608                         +          ( rkmjmim + rkmjpim +  rkmjmip + rkmjpip ) &
609                         + 4.0_wp * r(kp1,j,i)                                 &
610                         + 2.0_wp * ( r(kp1,j,i-1)   + r(kp1,j,i+1)   +        &
611                                      r(kp1,j+1,i)   + r(kp1,j-1,i)   )        &
612                         +          ( r(kp1,j-1,i-1) + r(kp1,j+1,i-1) +        &
613                                      r(kp1,j-1,i+1) + r(kp1,j+1,i+1) )        &
614                                                        )
615                ENDDO
616             ENDDO
617          ENDDO
618          !$OMP END PARALLEL
619
620       ENDIF
621
622!
623!--    Horizontal boundary conditions
624       CALL exchange_horiz( f_mg, 1)
625
626       IF ( .NOT. bc_lr_cyc )  THEN
627          IF (inflow_l .OR. outflow_l)  f_mg(:,:,nxl_mg(l)-1) = f_mg(:,:,nxl_mg(l))
628          IF (inflow_r .OR. outflow_r)  f_mg(:,:,nxr_mg(l)+1) = f_mg(:,:,nxr_mg(l))
629       ENDIF
630
631       IF ( .NOT. bc_ns_cyc )  THEN
632          IF (inflow_n .OR. outflow_n)  f_mg(:,nyn_mg(l)+1,:) = f_mg(:,nyn_mg(l),:)
633          IF (inflow_s .OR. outflow_s)  f_mg(:,nys_mg(l)-1,:) = f_mg(:,nys_mg(l),:)
634       ENDIF
635
636!
637!--    Boundary conditions at bottom and top of the domain.
638!--    These points are not handled by the above loop. Points may be within
639!--    buildings, but that doesn't matter.
640       IF ( ibc_p_b == 1 )  THEN
641          f_mg(nzb,:,: ) = f_mg(nzb+1,:,:)
642       ELSE
643          f_mg(nzb,:,: ) = 0.0_wp
644       ENDIF
645
646       IF ( ibc_p_t == 1 )  THEN
647          f_mg(nzt_mg(l)+1,:,: ) = f_mg(nzt_mg(l),:,:)
648       ELSE
649          f_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
650       ENDIF
651
652       CALL cpu_log( log_point_s(54), 'restrict', 'stop' )
653!
654!--    Why do we need a sorting here?
655       CALL sort_k_to_even_odd_blocks( f_mg , l)
656
657    END SUBROUTINE restrict_fast
658
659
660
661    SUBROUTINE prolong_fast( p, temp )
662
663!------------------------------------------------------------------------------!
664! Description:
665! ------------
666! Interpolates the correction of the perturbation pressure
667! to the next finer grid.
668!------------------------------------------------------------------------------!
669
670       USE control_parameters,                                                 &
671           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
672                  inflow_n, inflow_r, inflow_s, outflow_l, outflow_n,          &
673                  outflow_r, outflow_s
674
675       USE indices,                                                            &
676           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
677
678       IMPLICIT NONE
679
680       INTEGER(iwp) ::  i   !:
681       INTEGER(iwp) ::  j   !:
682       INTEGER(iwp) ::  k   !:
683       INTEGER(iwp) ::  l   !:
684       INTEGER(iwp) ::  kp1 !:
685       INTEGER(iwp) ::  ke  !: Index for prolog even
686       INTEGER(iwp) ::  ko  !: Index for prolog odd
687
688       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
689                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
690                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1 ) ::  p  !:
691
692       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
693                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
694                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  temp  !:
695
696
697       CALL cpu_log( log_point_s(55), 'prolong', 'start' )
698
699!
700!--    First, store elements of the coarser grid on the next finer grid
701       l = grid_level
702       ind_even_odd = even_odd_level(grid_level-1)
703
704       !$OMP PARALLEL PRIVATE (i,j,k,kp1,ke,ko)
705       !$OMP DO
706       DO  i = nxl_mg(l-1), nxr_mg(l-1)
707          DO  j = nys_mg(l-1), nyn_mg(l-1)
708
709             !DIR$ IVDEP
710             DO k = ind_even_odd+1, nzt_mg(l-1)
711                kp1 = k - ind_even_odd
712                ke  = 2 * ( k-ind_even_odd - 1 ) + 1
713                ko  = 2 * k - 1
714!
715!--             Points of the coarse grid are directly stored on the next finer
716!--             grid
717                temp(ko,2*j,2*i)   = p(k,j,i)
718!
719!--             Points between two coarse-grid points
720                temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
721                temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
722                temp(ke,2*j,2*i)   = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
723!
724!--             Points in the center of the planes stretched by four points
725!--             of the coarse grid cube
726                temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
727                                                   p(k,j+1,i) + p(k,j+1,i+1) )
728                temp(ke,2*j,2*i+1)   = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
729                                                   p(kp1,j,i) + p(kp1,j,i+1) )
730                temp(ke,2*j+1,2*i)   = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +   &
731                                                   p(kp1,j,i) + p(kp1,j+1,i) )
732!
733!--             Points in the middle of coarse grid cube
734                temp(ke,2*j+1,2*i+1) = 0.125_wp *                              &
735                                               ( p(k,j,i)     + p(k,j,i+1)   + &
736                                                 p(k,j+1,i)   + p(k,j+1,i+1) + &
737                                                 p(kp1,j,i)   + p(kp1,j,i+1) + &
738                                                 p(kp1,j+1,i) + p(kp1,j+1,i+1) )
739             ENDDO
740
741             !DIR$ IVDEP
742             DO k = nzb+1, ind_even_odd
743                kp1 = k + ind_even_odd + 1
744                ke  = 2 * k
745                ko  = 2 * ( k + ind_even_odd )
746!
747!--             Points of the coarse grid are directly stored on the next finer
748!--             grid
749                temp(ko,2*j,2*i)   = p(k,j,i)
750!
751!--             Points between two coarse-grid points
752                temp(ko,2*j,2*i+1) = 0.5_wp * ( p(k,j,i) + p(k,j,i+1) )
753                temp(ko,2*j+1,2*i) = 0.5_wp * ( p(k,j,i) + p(k,j+1,i) )
754                temp(ke,2*j,2*i)   = 0.5_wp * ( p(k,j,i) + p(kp1,j,i) )
755!
756!--             Points in the center of the planes stretched by four points
757!--             of the coarse grid cube
758                temp(ko,2*j+1,2*i+1) = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
759                                                   p(k,j+1,i) + p(k,j+1,i+1) )
760                temp(ke,2*j,2*i+1)   = 0.25_wp * ( p(k,j,i)   + p(k,j,i+1) +   &
761                                                   p(kp1,j,i) + p(kp1,j,i+1) )
762                temp(ke,2*j+1,2*i)   = 0.25_wp * ( p(k,j,i)   + p(k,j+1,i) +   &
763                                                   p(kp1,j,i) + p(kp1,j+1,i) )
764!
765!--             Points in the middle of coarse grid cube
766                temp(ke,2*j+1,2*i+1) = 0.125_wp *                              &
767                                               ( p(k,j,i)     + p(k,j,i+1)   + &
768                                                 p(k,j+1,i)   + p(k,j+1,i+1) + &
769                                                 p(kp1,j,i)   + p(kp1,j,i+1) + &
770                                                 p(kp1,j+1,i) + p(kp1,j+1,i+1) )
771             ENDDO
772
773          ENDDO
774       ENDDO
775       !$OMP END PARALLEL
776
777       ind_even_odd = even_odd_level(grid_level)
778!
779!--    Horizontal boundary conditions
780       CALL exchange_horiz( temp, 1)
781
782       IF ( .NOT. bc_lr_cyc )  THEN
783          IF (inflow_l .OR. outflow_l)  temp(:,:,nxl_mg(l)-1) = temp(:,:,nxl_mg(l))
784          IF (inflow_r .OR. outflow_r)  temp(:,:,nxr_mg(l)+1) = temp(:,:,nxr_mg(l))
785       ENDIF
786
787       IF ( .NOT. bc_ns_cyc )  THEN
788          IF (inflow_n .OR. outflow_n)  temp(:,nyn_mg(l)+1,:) = temp(:,nyn_mg(l),:)
789          IF (inflow_s .OR. outflow_s)  temp(:,nys_mg(l)-1,:) = temp(:,nys_mg(l),:)
790       ENDIF
791
792!
793!--    Bottom and top boundary conditions
794       IF ( ibc_p_b == 1 )  THEN
795          temp(nzb,:,: ) = temp(ind_even_odd+1,:,:)
796       ELSE
797          temp(nzb,:,: ) = 0.0_wp
798       ENDIF
799
800       IF ( ibc_p_t == 1 )  THEN
801          temp(nzt_mg(l)+1,:,: ) = temp(ind_even_odd,:,:)
802       ELSE
803          temp(nzt_mg(l)+1,:,: ) = 0.0_wp
804       ENDIF
805
806       CALL cpu_log( log_point_s(55), 'prolong', 'stop' )
807
808    END SUBROUTINE prolong_fast
809
810
811
812    SUBROUTINE redblack_fast( f_mg, p_mg )
813
814!------------------------------------------------------------------------------!
815! Description:
816! ------------
817! Relaxation method for the multigrid scheme. A Gauss-Seidel iteration with
818! 3D-Red-Black decomposition (GS-RB) is used.
819!------------------------------------------------------------------------------!
820
821       USE arrays_3d,                                                          &
822           ONLY:  f1_mg, f2_mg, f3_mg
823
824       USE control_parameters,                                                 &
825           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t, inflow_l,&
826                  inflow_n, inflow_r, inflow_s, masking_method, ngsrb,         &
827                  outflow_l, outflow_n, outflow_r, outflow_s, topography
828
829       USE grid_variables,                                                     &
830           ONLY:  ddx2_mg, ddy2_mg
831
832       USE indices,                                                            &
833           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
834                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
835                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
836                  nys_mg, nyn_mg, nzb, nzt_mg
837
838       IMPLICIT NONE
839
840       INTEGER(iwp) :: color    !:
841       INTEGER(iwp) :: i        !:
842       INTEGER(iwp) :: ic       !:
843       INTEGER(iwp) :: j        !:
844       INTEGER(iwp) :: jc       !:
845       INTEGER(iwp) :: jj       !:
846       INTEGER(iwp) :: k        !:
847       INTEGER(iwp) :: l        !:
848       INTEGER(iwp) :: n        !:
849       INTEGER(iwp) :: km1      !:
850       INTEGER(iwp) :: kp1      !:
851
852       LOGICAL :: unroll        !:
853
854       REAL(wp) ::  wall_left   !:
855       REAL(wp) ::  wall_north  !:
856       REAL(wp) ::  wall_right  !:
857       REAL(wp) ::  wall_south  !:
858       REAL(wp) ::  wall_total  !:
859       REAL(wp) ::  wall_top    !:
860
861       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
862                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
863                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f_mg  !:
864       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
865                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
866                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
867
868       l = grid_level
869
870!
871!--    Choose flag array of this level
872       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
873
874          SELECT CASE ( l )
875             CASE ( 1 )
876                flags => wall_flags_1
877             CASE ( 2 )
878                flags => wall_flags_2
879             CASE ( 3 )
880                flags => wall_flags_3
881             CASE ( 4 )
882                flags => wall_flags_4
883             CASE ( 5 )
884                flags => wall_flags_5
885             CASE ( 6 )
886                flags => wall_flags_6
887             CASE ( 7 )
888                flags => wall_flags_7
889             CASE ( 8 )
890                flags => wall_flags_8
891             CASE ( 9 )
892                flags => wall_flags_9
893             CASE ( 10 )
894                flags => wall_flags_10
895          END SELECT
896
897       ENDIF
898
899       unroll = ( MOD( nyn_mg(l)-nys_mg(l)+1, 4 ) == 0  .AND. &
900                  MOD( nxr_mg(l)-nxl_mg(l)+1, 2 ) == 0 )
901
902       DO  n = 1, ngsrb
903       
904          DO  color = 1, 2
905!
906!--          No wall treatment in case of flat topography
907             IF ( topography == 'flat'  .OR.  masking_method )  THEN
908
909                IF ( .NOT. unroll )  THEN
910
911                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'start' )
912!
913!--                Without unrolling of loops, no cache optimization
914                   !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
915                   !$OMP DO
916                   DO  i = nxl_mg(l), nxr_mg(l), 2
917                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
918                         !DIR$ IVDEP
919                         DO  k = ind_even_odd+1, nzt_mg(l)
920                            km1 = k-ind_even_odd-1
921                            kp1 = k-ind_even_odd
922                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
923                               ddx2_mg(l) *                                    &
924                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
925                               + ddy2_mg(l) *                                  &
926                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
927                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
928                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
929                               - f_mg(k,j,i)               )
930                         ENDDO
931                      ENDDO
932                   ENDDO
933   
934                   !$OMP DO
935                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
936                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
937                          !DIR$ IVDEP
938                          DO  k = ind_even_odd+1, nzt_mg(l)
939                            km1 = k-ind_even_odd-1
940                            kp1 = k-ind_even_odd
941                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
942                               ddx2_mg(l) *                                    &
943                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
944                               + ddy2_mg(l) *                                  &
945                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
946                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
947                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
948                               - f_mg(k,j,i)               )
949                         ENDDO
950                      ENDDO
951                   ENDDO
952 
953                   !$OMP DO
954                   DO  i = nxl_mg(l), nxr_mg(l), 2
955                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
956                         !DIR$ IVDEP
957                         DO  k = nzb+1, ind_even_odd
958                            km1 = k+ind_even_odd
959                            kp1 = k+ind_even_odd+1
960                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
961                               ddx2_mg(l) *                                    &
962                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
963                               + ddy2_mg(l) *                                  &
964                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
965                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
966                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
967                               - f_mg(k,j,i)               )
968                         ENDDO
969                      ENDDO
970                   ENDDO
971
972                   !$OMP DO
973                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
974                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
975                         !DIR$ IVDEP
976                         DO  k = nzb+1, ind_even_odd
977                            km1 = k+ind_even_odd
978                            kp1 = k+ind_even_odd+1
979                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
980                               ddx2_mg(l) *                                    &
981                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
982                               + ddy2_mg(l) *                                  &
983                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
984                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
985                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
986                               - f_mg(k,j,i)               )
987                         ENDDO
988                      ENDDO
989                   ENDDO
990                   !$OMP END PARALLEL
991
992                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll_f', 'stop' )
993
994                ELSE
995!
996!--                Loop unrolling along y, only one i loop for better cache use
997                   CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'start' )
998
999                   !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
1000                   !$OMP DO
1001                   DO  ic = nxl_mg(l), nxr_mg(l), 2
1002                      DO  jc = nys_mg(l), nyn_mg(l), 4
1003                         i  = ic
1004                         jj = jc+2-color
1005                         !DIR$ IVDEP
1006                         DO  k = ind_even_odd+1, nzt_mg(l)
1007                            km1 = k-ind_even_odd-1
1008                            kp1 = k-ind_even_odd
1009                            j   = jj
1010                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1011                               ddx2_mg(l) *                                    &
1012                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
1013                               + ddy2_mg(l) *                                  &
1014                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
1015                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1016                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1017                               - f_mg(k,j,i)               )
1018                            j = jj+2
1019                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1020                               ddx2_mg(l) *                                    &
1021                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
1022                               + ddy2_mg(l) *                                  &
1023                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
1024                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1025                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1026                               - f_mg(k,j,i)               )
1027                         ENDDO
1028
1029                         i  = ic+1
1030                         jj = jc+color-1
1031                         !DIR$ IVDEP
1032                         DO  k = ind_even_odd+1, nzt_mg(l)
1033                            km1 = k-ind_even_odd-1
1034                            kp1 = k-ind_even_odd
1035                            j   = jj
1036                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1037                               ddx2_mg(l) *                                    &
1038                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
1039                               + ddy2_mg(l) *                                  &
1040                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
1041                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1042                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1043                               - f_mg(k,j,i)               )
1044                            j = jj+2
1045                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1046                               ddx2_mg(l) *                                    &
1047                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
1048                               + ddy2_mg(l) *                                  &
1049                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
1050                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1051                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1052                               - f_mg(k,j,i)               )
1053                         ENDDO
1054
1055                         i  = ic
1056                         jj = jc+color-1
1057                         !DIR$ IVDEP
1058                         DO  k = nzb+1, ind_even_odd
1059                            km1 = k+ind_even_odd
1060                            kp1 = k+ind_even_odd+1
1061                            j   = jj
1062                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1063                               ddx2_mg(l) *                                    &
1064                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
1065                               + ddy2_mg(l) *                                  &
1066                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
1067                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1068                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1069                               - f_mg(k,j,i)               )
1070                            j = jj+2
1071                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1072                               ddx2_mg(l) *                                    &
1073                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
1074                               + ddy2_mg(l) *                                  &
1075                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
1076                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1077                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1078                               - f_mg(k,j,i)               )
1079                         ENDDO
1080
1081                         i  = ic+1
1082                         jj = jc+2-color
1083                         !DIR$ IVDEP
1084                         DO  k = nzb+1, ind_even_odd
1085                            km1 = k+ind_even_odd
1086                            kp1 = k+ind_even_odd+1
1087                            j   = jj
1088                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1089                               ddx2_mg(l) *                                    &
1090                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1)  )              &
1091                               + ddy2_mg(l) *                                  &
1092                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i) )               &
1093                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1094                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1095                               - f_mg(k,j,i)               )
1096                            j = jj+2
1097                            p_mg(k,j,i) = 1.0_wp / f1_mg_b(k,l) * (            &
1098                               ddx2_mg(l) *                                    &
1099                               ( p_mg(k,j,i+1) + p_mg(k,j,i-1) )               &
1100                               + ddy2_mg(l) *                                  &
1101                               ( p_mg(k,j+1,i) + p_mg(k,j-1,i)  )              &
1102                               + f2_mg_b(k,l) * p_mg(kp1,j,i)                  &
1103                               + f3_mg_b(k,l) * p_mg(km1,j,i)                  &
1104                               - f_mg(k,j,i)               )
1105                         ENDDO
1106
1107                      ENDDO
1108                   ENDDO
1109                   !$OMP END PARALLEL
1110
1111                   CALL cpu_log( log_point_s(38), 'redblack_unroll_f', 'stop' )
1112
1113                ENDIF
1114
1115             ELSE
1116
1117                IF ( .NOT. unroll )  THEN
1118
1119                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'start' )
1120!
1121!--                Without unrolling of loops, no cache optimization /
1122!--                vectorization. Therefore, the non-vectorized IBITS function
1123!--                is still used.
1124                   !$OMP PARALLEL PRIVATE (i,j,k,km1,kp1)
1125                   !$OMP DO
1126                   DO  i = nxl_mg(l), nxr_mg(l), 2
1127                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1128                         !DIR$ IVDEP
1129                         DO  k = ind_even_odd+1, nzt_mg(l)
1130                            km1 = k-ind_even_odd-1
1131                            kp1 = k-ind_even_odd
1132
1133                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
1134                               ddx2_mg(l) *                                    &
1135                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1136                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
1137                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
1138                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
1139                               + ddy2_mg(l) *                                  &
1140                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1141                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
1142                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
1143                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
1144                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
1145                               + f3_mg(k,l) *                                  &
1146                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1147                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
1148                               - f_mg(k,j,i)               )
1149                         ENDDO
1150                      ENDDO
1151                   ENDDO
1152   
1153                   !$OMP DO
1154                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
1155                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
1156                          !DIR$ IVDEP
1157                          DO  k = ind_even_odd+1, nzt_mg(l)
1158                            km1 = k-ind_even_odd-1
1159                            kp1 = k-ind_even_odd
1160
1161                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
1162                               ddx2_mg(l) *                                    &
1163                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1164                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
1165                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
1166                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
1167                               + ddy2_mg(l) *                                  &
1168                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1169                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
1170                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
1171                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
1172                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
1173                               + f3_mg(k,l) *                                  &
1174                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1175                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
1176                               - f_mg(k,j,i)               )
1177                         ENDDO
1178                      ENDDO
1179                   ENDDO
1180
1181                   !$OMP DO
1182                   DO  i = nxl_mg(l), nxr_mg(l), 2
1183                      DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
1184                         !DIR$ IVDEP
1185                         DO  k = nzb+1, ind_even_odd
1186                            km1 = k+ind_even_odd
1187                            kp1 = k+ind_even_odd+1
1188                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
1189                               ddx2_mg(l) *                                    &
1190                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1191                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
1192                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
1193                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
1194                               + ddy2_mg(l) *                                  &
1195                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1196                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
1197                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
1198                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
1199                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
1200                               + f3_mg(k,l) *                                  &
1201                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1202                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) ) &
1203                               - f_mg(k,j,i)               )
1204                         ENDDO
1205                      ENDDO
1206                   ENDDO
1207
1208                   !$OMP DO
1209                   DO  i = nxl_mg(l)+1, nxr_mg(l), 2
1210                      DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
1211                         !DIR$ IVDEP
1212                         DO  k = nzb+1, ind_even_odd
1213                            km1 = k+ind_even_odd
1214                            kp1 = k+ind_even_odd+1
1215
1216                            p_mg(k,j,i) = 1.0_wp / f1_mg(k,l) * (              &
1217                               ddx2_mg(l) *                                    &
1218                               ( p_mg(k,j,i+1) + IBITS( flags(k,j,i), 5, 1 ) * &
1219                               ( p_mg(k,j,i) - p_mg(k,j,i+1) ) +               &
1220                               p_mg(k,j,i-1) + IBITS( flags(k,j,i), 4, 1 ) *   &
1221                               ( p_mg(k,j,i) - p_mg(k,j,i-1) ) )               &
1222                               + ddy2_mg(l) *                                  &
1223                               ( p_mg(k,j+1,i) + IBITS( flags(k,j,i), 3, 1 ) * &
1224                               ( p_mg(k,j,i) - p_mg(k,j+1,i) ) +               &
1225                               p_mg(k,j-1,i) + IBITS( flags(k,j,i), 2, 1 ) *   &
1226                               ( p_mg(k,j,i) - p_mg(k,j-1,i) ) )               &
1227                               + f2_mg(k,l) * p_mg(kp1,j,i)                    &
1228                               + f3_mg(k,l) *                                  &
1229                               ( p_mg(km1,j,i) + IBITS( flags(k,j,i), 0, 1 ) * &
1230                               ( p_mg(k,j,i) - p_mg(km1,j,i) ) )               &
1231                               - f_mg(k,j,i)               )
1232                         ENDDO
1233                      ENDDO
1234                   ENDDO
1235                   !$OMP END PARALLEL
1236
1237                   CALL cpu_log( log_point_s(36), 'redblack_no_unroll', 'stop' )
1238
1239                ELSE
1240!
1241!--                Loop unrolling along y, only one i loop for better cache use
1242                   CALL cpu_log( log_point_s(38), 'redblack_unroll', 'start' )
1243
1244                   !$OMP PARALLEL PRIVATE (i,j,k,ic,jc,km1,kp1,jj)
1245                   !$OMP DO
1246                   DO  ic = nxl_mg(l), nxr_mg(l), 2
1247                      DO  jc = nys_mg(l), nyn_mg(l), 4
1248                         i  = ic
1249                         jj = jc+2-color
1250                         !DIR$ IVDEP
1251                         DO  k = ind_even_odd+1, nzt_mg(l)
1252
1253                            km1 = k-ind_even_odd-1
1254                            kp1 = k-ind_even_odd
1255
1256                            j = jj
1257                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1258                                ddx2_mg(l) *                                                 &
1259                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1260                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1261                                + ddy2_mg(l) *                                               &
1262                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1263                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1264                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1265                                + f3_mg(k,l) *                                               &
1266                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1267                                - f_mg(k,j,i)               )
1268                            j = jj+2
1269                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1270                                ddx2_mg(l) *                                                 &
1271                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1272                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1273                                + ddy2_mg(l) *                                               &
1274                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1275                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1276                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1277                                + f3_mg(k,l) *                                               &
1278                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1279                                - f_mg(k,j,i)               )
1280
1281                         ENDDO
1282
1283                         i  = ic+1
1284                         jj = jc+color-1
1285                         !DIR$ IVDEP
1286                         DO  k = ind_even_odd+1, nzt_mg(l)
1287
1288                            km1 = k-ind_even_odd-1
1289                            kp1 = k-ind_even_odd
1290
1291                            j =jj
1292                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1293                                ddx2_mg(l) *                                                 &
1294                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1295                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1296                                + ddy2_mg(l) *                                               &
1297                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1298                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1299                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1300                                + f3_mg(k,l) *                                               &
1301                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1302                                - f_mg(k,j,i)               )
1303
1304                            j = jj+2
1305                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1306                                ddx2_mg(l) *                                                 &
1307                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1308                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1309                                + ddy2_mg(l) *                                               &
1310                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1311                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1312                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1313                                + f3_mg(k,l) *                                               &
1314                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1315                                - f_mg(k,j,i)               )
1316
1317                         ENDDO
1318
1319                         i  = ic
1320                         jj = jc+color-1
1321                         !DIR$ IVDEP
1322                         DO  k = nzb+1, ind_even_odd
1323
1324                            km1 = k+ind_even_odd
1325                            kp1 = k+ind_even_odd+1
1326
1327                            j =jj
1328                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1329                                ddx2_mg(l) *                                                 &
1330                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1331                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1332                                + ddy2_mg(l) *                                               &
1333                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1334                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1335                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1336                                + f3_mg(k,l) *                                               &
1337                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1338                                - f_mg(k,j,i)               )
1339
1340                            j = jj+2
1341                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1342                                ddx2_mg(l) *                                                 &
1343                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1344                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1345                                + ddy2_mg(l) *                                               &
1346                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1347                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1348                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1349                                + f3_mg(k,l) *                                               &
1350                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1351                                - f_mg(k,j,i)               )
1352
1353                         ENDDO
1354
1355                         i  = ic+1
1356                         jj = jc+2-color
1357                         !DIR$ IVDEP
1358                         DO  k = nzb+1, ind_even_odd
1359
1360                            km1 = k+ind_even_odd
1361                            kp1 = k+ind_even_odd+1
1362
1363                            j =jj
1364                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1365                                ddx2_mg(l) *                                                 &
1366                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1367                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1368                                + ddy2_mg(l) *                                               &
1369                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1370                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1371                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1372                                + f3_mg(k,l) *                                               &
1373                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1374                                - f_mg(k,j,i)               )
1375
1376                            j = jj+2
1377                            p_mg(k,j,i) = 1.0 / f1_mg(k,l) * (                               &
1378                                ddx2_mg(l) *                                                 &
1379                                ( MERGE(p_mg(k,j,i),p_mg(k,j,i+1),BTEST( flags(k,j,i), 5))   &
1380                                + MERGE(p_mg(k,j,i),p_mg(k,j,i-1),BTEST( flags(k,j,i), 4)) ) &
1381                                + ddy2_mg(l) *                                               &
1382                                ( MERGE(p_mg(k,j,i),p_mg(k,j+1,i),BTEST( flags(k,j,i), 3))   &
1383                                + MERGE(p_mg(k,j,i),p_mg(k,j-1,i),BTEST( flags(k,j,i), 2)) ) &
1384                                + f2_mg(k,l) * p_mg(kp1,j,i)                                 &
1385                                + f3_mg(k,l) *                                               &
1386                                MERGE(p_mg(k,j,i),p_mg(km1,j,i),BTEST( flags(k,j,i), 0))     &
1387                                - f_mg(k,j,i)               )
1388
1389                         ENDDO
1390
1391                      ENDDO
1392                   ENDDO
1393                   !$OMP END PARALLEL
1394
1395                   CALL cpu_log( log_point_s(38), 'redblack_unroll', 'stop' )
1396
1397                ENDIF
1398
1399             ENDIF
1400
1401!
1402!--          Horizontal boundary conditions
1403             CALL special_exchange_horiz( p_mg, color )
1404
1405              IF ( .NOT. bc_lr_cyc )  THEN
1406                IF ( inflow_l .OR. outflow_l )  THEN
1407                   p_mg(:,:,nxl_mg(l)-1) = p_mg(:,:,nxl_mg(l))
1408                ENDIF
1409                IF ( inflow_r .OR. outflow_r )  THEN
1410                   p_mg(:,:,nxr_mg(l)+1) = p_mg(:,:,nxr_mg(l))
1411                ENDIF
1412             ENDIF
1413
1414             IF ( .NOT. bc_ns_cyc )  THEN
1415                IF ( inflow_n .OR. outflow_n )  THEN
1416                   p_mg(:,nyn_mg(l)+1,:) = p_mg(:,nyn_mg(l),:)
1417                ENDIF
1418                IF ( inflow_s .OR. outflow_s )  THEN
1419                   p_mg(:,nys_mg(l)-1,:) = p_mg(:,nys_mg(l),:)
1420                ENDIF
1421             ENDIF
1422
1423!
1424!--          Bottom and top boundary conditions
1425             IF ( ibc_p_b == 1 )  THEN
1426                p_mg(nzb,:,: ) = p_mg(ind_even_odd+1,:,:)
1427             ELSE
1428                p_mg(nzb,:,: ) = 0.0_wp
1429             ENDIF
1430
1431             IF ( ibc_p_t == 1 )  THEN
1432                p_mg(nzt_mg(l)+1,:,: ) = p_mg(ind_even_odd,:,:)
1433             ELSE
1434                p_mg(nzt_mg(l)+1,:,: ) = 0.0_wp
1435             ENDIF
1436
1437          ENDDO
1438       ENDDO
1439!
1440!--    Set pressure within topography and at the topography surfaces
1441       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
1442
1443          !$OMP PARALLEL PRIVATE (i,j,k,kp1,wall_left,wall_north,wall_right,   &
1444          !$OMP                   wall_south,wall_top,wall_total)
1445          !$OMP DO
1446          DO  i = nxl_mg(l), nxr_mg(l)
1447             DO  j = nys_mg(l), nyn_mg(l)
1448                !DIR$ IVDEP
1449                DO  k = ind_even_odd+1, nzt_mg(l)
1450                   km1 = k-ind_even_odd-1
1451                   kp1 = k-ind_even_odd
1452!
1453!--                First, set pressure inside topography to zero
1454                   p_mg(k,j,i) = p_mg(k,j,i) *                                 &
1455                                 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1456!
1457!--                Second, determine if the gridpoint inside topography is
1458!--                adjacent to a wall and set its value to a value given by the
1459!--                average of those values obtained from Neumann boundary
1460!--                condition
1461                   wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
1462                   wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1463                   wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1464                   wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1465                   wall_top   = IBITS( flags(kp1,j,i), 0, 1 )
1466                   wall_total = wall_left + wall_right + wall_south +          &
1467                                wall_north + wall_top
1468
1469                   IF ( wall_total > 0.0_wp )  THEN
1470                      p_mg(k,j,i) = 1.0_wp / wall_total *                      &
1471                                    ( wall_left  * p_mg(k,j,i-1) +             &
1472                                      wall_right * p_mg(k,j,i+1) +             &
1473                                      wall_south * p_mg(k,j-1,i) +             &
1474                                      wall_north * p_mg(k,j+1,i) +             &
1475                                      wall_top   * p_mg(kp1,j,i) )
1476                   ENDIF
1477                ENDDO
1478
1479                !DIR$ IVDEP
1480                DO  k = nzb+1, ind_even_odd
1481                   km1 = k+ind_even_odd
1482                   kp1 = k+ind_even_odd+1
1483!
1484!--                First, set pressure inside topography to zero
1485                   p_mg(k,j,i) = p_mg(k,j,i) *                                 &
1486                                 ( 1.0_wp - IBITS( flags(k,j,i), 6, 1 ) )
1487!
1488!--                Second, determine if the gridpoint inside topography is
1489!--                adjacent to a wall and set its value to a value given by the
1490!--                average of those values obtained from Neumann boundary
1491!--                condition
1492                   wall_left  = IBITS( flags(k,j,i-1), 5, 1 )
1493                   wall_right = IBITS( flags(k,j,i+1), 4, 1 )
1494                   wall_south = IBITS( flags(k,j-1,i), 3, 1 )
1495                   wall_north = IBITS( flags(k,j+1,i), 2, 1 )
1496                   wall_top   = IBITS( flags(kp1,j,i), 0, 1 )
1497                   wall_total = wall_left + wall_right + wall_south +          &
1498                                wall_north + wall_top
1499
1500                   IF ( wall_total > 0.0_wp )  THEN
1501                      p_mg(k,j,i) = 1.0_wp / wall_total *                      &
1502                                    ( wall_left  * p_mg(k,j,i-1) +             &
1503                                      wall_right * p_mg(k,j,i+1) +             &
1504                                      wall_south * p_mg(k,j-1,i) +             &
1505                                      wall_north * p_mg(k,j+1,i) +             &
1506                                      wall_top   * p_mg(kp1,j,i) )
1507                   ENDIF
1508                ENDDO
1509             ENDDO
1510          ENDDO
1511          !$OMP END PARALLEL
1512
1513!
1514!--       One more time horizontal boundary conditions
1515          CALL exchange_horiz( p_mg, 1)
1516
1517       ENDIF
1518
1519    END SUBROUTINE redblack_fast
1520
1521
1522
1523    SUBROUTINE sort_k_to_even_odd_blocks( p_mg , glevel )
1524
1525!------------------------------------------------------------------------------!
1526! Description:
1527! ------------
1528! Sort k-Dimension from sequential into blocks of even and odd.
1529! This is required to vectorize the red-black subroutine.
1530! Version for 3D-REAL arrays
1531!------------------------------------------------------------------------------!
1532
1533       USE control_parameters,                                                 &
1534           ONLY:  grid_level
1535
1536       USE indices,                                                            &
1537           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1538
1539       IMPLICIT NONE
1540
1541       INTEGER(iwp), INTENT(IN) ::  glevel
1542
1543       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1,                               &
1544                           nys_mg(glevel)-1:nyn_mg(glevel)+1,                  &
1545                           nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::  p_mg  !:
1546!
1547!--    Local variables
1548       INTEGER(iwp) :: i        !:
1549       INTEGER(iwp) :: j        !:
1550       INTEGER(iwp) :: k        !:
1551       INTEGER(iwp) :: l        !:
1552       INTEGER(iwp) :: ind      !:
1553       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp  !:
1554
1555
1556       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
1557
1558       l = glevel
1559       ind_even_odd = even_odd_level(l)
1560
1561       !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
1562       !$OMP DO
1563       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
1564          DO  j = nys_mg(l)-1, nyn_mg(l)+1
1565
1566!
1567!--          Sort the data with even k index
1568             ind = nzb-1
1569             DO  k = nzb, nzt_mg(l), 2
1570                ind = ind + 1
1571                tmp(ind) = p_mg(k,j,i)
1572             ENDDO
1573!
1574!--          Sort the data with odd k index
1575             DO  k = nzb+1, nzt_mg(l)+1, 2
1576                ind = ind + 1
1577                tmp(ind) = p_mg(k,j,i)
1578             ENDDO
1579
1580             p_mg(:,j,i) = tmp
1581
1582          ENDDO
1583       ENDDO
1584       !$OMP END PARALLEL
1585
1586       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
1587
1588    END SUBROUTINE sort_k_to_even_odd_blocks
1589
1590
1591
1592    SUBROUTINE sort_k_to_even_odd_blocks_1d( f_mg, f_mg_b, glevel )
1593
1594!------------------------------------------------------------------------------!
1595! Description:
1596! ------------
1597! Sort k-Dimension from sequential into blocks of even and odd.
1598! This is required to vectorize the red-black subroutine.
1599! Version for 1D-REAL arrays
1600!------------------------------------------------------------------------------!
1601
1602       USE indices,                                                            &
1603           ONLY:  nzb, nzt_mg
1604
1605       IMPLICIT NONE
1606
1607       INTEGER(iwp), INTENT(IN) ::  glevel
1608
1609       REAL(wp), DIMENSION(nzb+1:nzt_mg(glevel)) ::  f_mg
1610       REAL(wp), DIMENSION(nzb:nzt_mg(glevel)+1) ::  f_mg_b
1611
1612!
1613!--    Local variables
1614       INTEGER(iwp) :: ind   !:
1615       INTEGER(iwp) :: k     !:
1616
1617
1618       ind = nzb - 1
1619!
1620!--    Sort the data with even k index
1621       DO  k = nzb, nzt_mg(glevel), 2
1622          ind = ind + 1
1623          IF ( k >= nzb+1  .AND.  k <= nzt_mg(glevel) )  THEN
1624             f_mg_b(ind) = f_mg(k)
1625          ENDIF
1626       ENDDO
1627!
1628!--    Sort the data with odd k index
1629       DO  k = nzb+1, nzt_mg(glevel)+1, 2
1630          ind = ind + 1
1631          IF( k >= nzb+1  .AND.  k <= nzt_mg(glevel) )  THEN
1632             f_mg_b(ind) = f_mg(k)
1633          ENDIF
1634       ENDDO
1635
1636    END SUBROUTINE sort_k_to_even_odd_blocks_1d
1637
1638
1639
1640    SUBROUTINE sort_k_to_even_odd_blocks_int( i_mg , glevel )
1641
1642!------------------------------------------------------------------------------!
1643! Description:
1644! ------------
1645! Sort k-Dimension from sequential into blocks of even and odd.
1646! This is required to vectorize the red-black subroutine.
1647! Version for 2D-INTEGER arrays
1648!------------------------------------------------------------------------------!
1649
1650       USE control_parameters,                                                 &
1651           ONLY:  grid_level
1652
1653       USE indices,                                                            &
1654           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1655
1656       IMPLICIT NONE
1657
1658       INTEGER(iwp), INTENT(IN) ::  glevel
1659
1660       INTEGER(iwp), DIMENSION(nzb:nzt_mg(glevel)+1,                           &
1661                               nys_mg(glevel)-1:nyn_mg(glevel)+1,              &
1662                               nxl_mg(glevel)-1:nxr_mg(glevel)+1) ::  i_mg  !:
1663!
1664!--    Local variables
1665       INTEGER(iwp) :: i        !:
1666       INTEGER(iwp) :: j        !:
1667       INTEGER(iwp) :: k        !:
1668       INTEGER(iwp) :: l        !:
1669       INTEGER(iwp) :: ind      !:
1670       INTEGER(iwp),DIMENSION(nzb:nzt_mg(glevel)+1) ::  tmp
1671
1672
1673       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'start' )
1674
1675       l = glevel
1676       ind_even_odd = even_odd_level(l)
1677
1678       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
1679          DO  j = nys_mg(l)-1, nyn_mg(l)+1
1680
1681!
1682!--          Sort the data with even k index
1683             ind = nzb-1
1684             DO  k = nzb, nzt_mg(l), 2
1685                ind = ind + 1
1686                tmp(ind) = i_mg(k,j,i)
1687             ENDDO
1688!
1689!++          ATTENTION: Check reason for this error. Remove it or replace WRITE
1690!++                     by PALM message
1691#if defined ( __parallel )
1692             IF ( ind /= ind_even_odd )  THEN
1693                WRITE (0,*) 'ERROR ==> illegal ind_even_odd ',ind,ind_even_odd,l
1694                CALL MPI_ABORT(MPI_COMM_WORLD,i,j)
1695             ENDIF
1696#endif
1697!
1698!--          Sort the data with odd k index
1699             DO  k = nzb+1, nzt_mg(l)+1, 2
1700                ind = ind + 1
1701                tmp(ind) = i_mg(k,j,i)
1702             ENDDO
1703
1704             i_mg(:,j,i) = tmp
1705
1706          ENDDO
1707       ENDDO
1708
1709       CALL cpu_log( log_point_s(52), 'sort_k_to_even_odd', 'stop' )
1710
1711    END SUBROUTINE sort_k_to_even_odd_blocks_int
1712
1713
1714
1715    SUBROUTINE sort_k_to_sequential( p_mg )
1716
1717!------------------------------------------------------------------------------!
1718! Description:
1719! ------------
1720! Sort k-dimension from blocks of even and odd into sequential
1721!------------------------------------------------------------------------------!
1722
1723       USE control_parameters,                                                 &
1724           ONLY:  grid_level
1725
1726       USE indices,                                                            &
1727           ONLY:  nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1728
1729       IMPLICIT NONE
1730
1731       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1732                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1733                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
1734!
1735!--    Local variables
1736       INTEGER(iwp) :: i        !:
1737       INTEGER(iwp) :: j        !:
1738       INTEGER(iwp) :: k        !:
1739       INTEGER(iwp) :: l        !:
1740       INTEGER(iwp) :: ind      !:
1741
1742       REAL(wp),DIMENSION(nzb:nzt_mg(grid_level)+1) ::  tmp
1743
1744
1745       l = grid_level
1746
1747       !$OMP PARALLEL PRIVATE (i,j,k,ind,tmp)
1748       !$OMP DO
1749       DO  i = nxl_mg(l)-1, nxr_mg(l)+1
1750          DO  j = nys_mg(l)-1, nyn_mg(l)+1
1751
1752             ind = nzb - 1
1753             tmp = p_mg(:,j,i)
1754             DO  k = nzb, nzt_mg(l), 2
1755                ind = ind + 1
1756                p_mg(k,j,i) = tmp(ind)
1757             ENDDO
1758
1759             DO  k = nzb+1, nzt_mg(l)+1, 2
1760                ind = ind + 1
1761                p_mg(k,j,i) = tmp(ind)
1762             ENDDO
1763          ENDDO
1764       ENDDO
1765       !$OMP END PARALLEL
1766
1767    END SUBROUTINE sort_k_to_sequential
1768
1769
1770
1771    SUBROUTINE mg_gather_fast( f2, f2_sub )
1772
1773       USE control_parameters,                                                 &
1774           ONLY:  grid_level
1775
1776       USE cpulog,                                                             &
1777           ONLY:  cpu_log, log_point_s
1778
1779       USE indices,                                                            &
1780           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1781
1782       IMPLICIT NONE
1783
1784       INTEGER(iwp) ::  i       !:
1785       INTEGER(iwp) ::  il      !:
1786       INTEGER(iwp) ::  ir      !:
1787       INTEGER(iwp) ::  j       !:
1788       INTEGER(iwp) ::  jn      !:
1789       INTEGER(iwp) ::  js      !:
1790       INTEGER(iwp) ::  k       !:
1791       INTEGER(iwp) ::  nwords  !:
1792
1793       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1794                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
1795                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2    !:
1796       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1797                       nys_mg(grid_level)-1:nyn_mg(grid_level)+1,              &
1798                       nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  f2_l  !:
1799
1800       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
1801                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
1802                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  f2_sub  !:
1803
1804
1805#if defined( __parallel )
1806       CALL cpu_log( log_point_s(34), 'mg_gather', 'start' )
1807
1808       f2_l = 0.0_wp
1809
1810!
1811!--    Store the local subdomain array on the total array
1812       js = mg_loc_ind(3,myid)
1813       IF ( south_border_pe )  js = js - 1
1814       jn = mg_loc_ind(4,myid)
1815       IF ( north_border_pe )  jn = jn + 1
1816       il = mg_loc_ind(1,myid)
1817       IF ( left_border_pe )   il = il - 1
1818       ir = mg_loc_ind(2,myid)
1819       IF ( right_border_pe )  ir = ir + 1
1820       DO  i = il, ir
1821          DO  j = js, jn
1822             DO  k = nzb, nzt_mg(grid_level)+1
1823                f2_l(k,j,i) = f2_sub(k,j,i)
1824             ENDDO
1825          ENDDO
1826       ENDDO
1827
1828!
1829!--    Find out the number of array elements of the total array
1830       nwords = SIZE( f2 )
1831
1832!
1833!--    Gather subdomain data from all PEs
1834       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1835       CALL MPI_ALLREDUCE( f2_l(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1), &
1836                           f2(nzb,nys_mg(grid_level)-1,nxl_mg(grid_level)-1),   &
1837                           nwords, MPI_REAL, MPI_SUM, comm2d, ierr )
1838
1839       CALL cpu_log( log_point_s(34), 'mg_gather', 'stop' )
1840#endif
1841   
1842    END SUBROUTINE mg_gather_fast
1843
1844
1845
1846    SUBROUTINE mg_scatter_fast( p2, p2_sub )
1847!
1848!--    TODO: It might be possible to improve the speed of this routine by using
1849!--    non-blocking communication
1850
1851       USE control_parameters,                                                 &
1852           ONLY:  grid_level
1853
1854       USE cpulog,                                                             &
1855           ONLY:  cpu_log, log_point_s
1856
1857       USE indices,                                                            &
1858           ONLY:  mg_loc_ind, nxl_mg, nxr_mg, nys_mg, nyn_mg, nzb, nzt_mg
1859
1860       IMPLICIT NONE
1861
1862       INTEGER(iwp) ::  nwords  !:
1863
1864       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1865                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
1866                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !:
1867
1868       REAL(wp), DIMENSION(nzb:mg_loc_ind(5,myid)+1,                           &
1869                           mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,          &
1870                           mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) ::  p2_sub  !:
1871
1872!
1873!--    Find out the number of array elements of the subdomain array
1874       nwords = SIZE( p2_sub )
1875
1876#if defined( __parallel )
1877       CALL cpu_log( log_point_s(35), 'mg_scatter', 'start' )
1878
1879       p2_sub = p2(:,mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1, &
1880                     mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1)
1881
1882       CALL cpu_log( log_point_s(35), 'mg_scatter', 'stop' )
1883#endif
1884   
1885    END SUBROUTINE mg_scatter_fast
1886
1887
1888
1889    RECURSIVE SUBROUTINE next_mg_level_fast( f_mg, p_mg, p3, r )
1890
1891!------------------------------------------------------------------------------!
1892! Description:
1893! ------------
1894! This is where the multigrid technique takes place. V- and W- Cycle are
1895! implemented and steered by the parameter "gamma". Parameter "nue" determines
1896! the convergence of the multigrid iterative solution. There are nue times
1897! RB-GS iterations. It should be set to "1" or "2", considering the time effort
1898! one would like to invest. Last choice shows a very good converging factor,
1899! but leads to an increase in computing time.
1900!------------------------------------------------------------------------------!
1901
1902       USE control_parameters,                                                 &
1903           ONLY:  bc_lr_dirrad, bc_lr_raddir, bc_ns_dirrad, bc_ns_raddir,      &
1904                  gamma_mg, grid_level, grid_level_count, ibc_p_b, ibc_p_t,    &
1905                  inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level,  &
1906                  mg_switch_to_pe0_level, mg_switch_to_pe0, ngsrb, outflow_l,  &
1907                  outflow_n, outflow_r, outflow_s
1908
1909       USE indices,                                                            &
1910           ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,      &
1911                  nyn_mg, nzb, nzt, nzt_mg
1912
1913       IMPLICIT NONE
1914
1915       INTEGER(iwp) ::  i            !:
1916       INTEGER(iwp) ::  j            !:
1917       INTEGER(iwp) ::  k            !:
1918       INTEGER(iwp) ::  nxl_mg_save  !:
1919       INTEGER(iwp) ::  nxr_mg_save  !:
1920       INTEGER(iwp) ::  nyn_mg_save  !:
1921       INTEGER(iwp) ::  nys_mg_save  !:
1922       INTEGER(iwp) ::  nzt_mg_save  !:
1923
1924       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1925                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1926                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: f_mg  !:
1927       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1928                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1929                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p_mg  !:
1930       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1931                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1932                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: p3    !:
1933       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
1934                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
1935                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) :: r     !:
1936
1937       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1938                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
1939                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  f2  !:
1940       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
1941                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
1942                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  p2  !:
1943
1944       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  f2_sub  !:
1945       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  p2_sub  !:
1946
1947!
1948!--    Restriction to the coarsest grid
1949    10 IF ( grid_level == 1 )  THEN
1950
1951!
1952!--       Solution on the coarsest grid. Double the number of Gauss-Seidel
1953!--       iterations in order to get a more accurate solution.
1954          ngsrb = 2 * ngsrb
1955
1956          ind_even_odd = even_odd_level(grid_level)
1957
1958          CALL redblack_fast( f_mg, p_mg )
1959
1960          ngsrb = ngsrb / 2
1961
1962
1963       ELSEIF ( grid_level /= 1 )  THEN
1964
1965          grid_level_count(grid_level) = grid_level_count(grid_level) + 1
1966
1967!
1968!--       Solution on the actual grid level
1969          ind_even_odd = even_odd_level(grid_level)
1970
1971          CALL redblack_fast( f_mg, p_mg )
1972
1973!
1974!--       Determination of the actual residual
1975          CALL resid_fast( f_mg, p_mg, r )
1976
1977!--       Restriction of the residual (finer grid values!) to the next coarser
1978!--       grid. Therefore, the grid level has to be decremented now. nxl..nzt have
1979!--       to be set to the coarse grid values, because these variables are needed
1980!--       for the exchange of ghost points in routine exchange_horiz
1981          grid_level = grid_level - 1
1982
1983          nxl = nxl_mg(grid_level)
1984          nys = nys_mg(grid_level)
1985          nxr = nxr_mg(grid_level)
1986          nyn = nyn_mg(grid_level)
1987          nzt = nzt_mg(grid_level)
1988
1989          IF ( grid_level == mg_switch_to_pe0_level )  THEN
1990
1991!
1992!--          From this level on, calculations are done on PE0 only.
1993!--          First, carry out restriction on the subdomain.
1994!--          Therefore, indices of the level have to be changed to subdomain
1995!--          values in between (otherwise, the restrict routine would expect
1996!--          the gathered array)
1997
1998             nxl_mg_save = nxl_mg(grid_level)
1999             nxr_mg_save = nxr_mg(grid_level)
2000             nys_mg_save = nys_mg(grid_level)
2001             nyn_mg_save = nyn_mg(grid_level)
2002             nzt_mg_save = nzt_mg(grid_level)
2003             nxl_mg(grid_level) = mg_loc_ind(1,myid)
2004             nxr_mg(grid_level) = mg_loc_ind(2,myid)
2005             nys_mg(grid_level) = mg_loc_ind(3,myid)
2006             nyn_mg(grid_level) = mg_loc_ind(4,myid)
2007             nzt_mg(grid_level) = mg_loc_ind(5,myid)
2008             nxl = mg_loc_ind(1,myid)
2009             nxr = mg_loc_ind(2,myid)
2010             nys = mg_loc_ind(3,myid)
2011             nyn = mg_loc_ind(4,myid)
2012             nzt = mg_loc_ind(5,myid)
2013
2014             ALLOCATE( f2_sub(nzb:nzt_mg(grid_level)+1,                    &
2015                              nys_mg(grid_level)-1:nyn_mg(grid_level)+1,   &
2016                              nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) )
2017
2018             CALL restrict_fast( f2_sub, r )
2019
2020!
2021!--          Restore the correct indices of this level
2022             nxl_mg(grid_level) = nxl_mg_save
2023             nxr_mg(grid_level) = nxr_mg_save
2024             nys_mg(grid_level) = nys_mg_save
2025             nyn_mg(grid_level) = nyn_mg_save
2026             nzt_mg(grid_level) = nzt_mg_save
2027             nxl = nxl_mg(grid_level)
2028             nxr = nxr_mg(grid_level)
2029             nys = nys_mg(grid_level)
2030             nyn = nyn_mg(grid_level)
2031             nzt = nzt_mg(grid_level)
2032!
2033!--          Gather all arrays from the subdomains on PE0
2034             CALL mg_gather_fast( f2, f2_sub )
2035
2036!
2037!--          Set switch for routine exchange_horiz, that no ghostpoint exchange
2038!--          has to be carried out from now on
2039             mg_switch_to_pe0 = .TRUE.
2040
2041!
2042!--          In case of non-cyclic lateral boundary conditions, both in- and
2043!--          outflow conditions have to be used on all PEs after the switch,
2044!--          because then they have the total domain.
2045             IF ( bc_lr_dirrad )  THEN
2046                inflow_l  = .TRUE.
2047                inflow_r  = .FALSE.
2048                outflow_l = .FALSE.
2049                outflow_r = .TRUE.
2050             ELSEIF ( bc_lr_raddir )  THEN
2051                inflow_l  = .FALSE.
2052                inflow_r  = .TRUE.
2053                outflow_l = .TRUE.
2054                outflow_r = .FALSE.
2055             ENDIF
2056
2057             IF ( bc_ns_dirrad )  THEN
2058                inflow_n  = .TRUE.
2059                inflow_s  = .FALSE.
2060                outflow_n = .FALSE.
2061                outflow_s = .TRUE.
2062             ELSEIF ( bc_ns_raddir )  THEN
2063                inflow_n  = .FALSE.
2064                inflow_s  = .TRUE.
2065                outflow_n = .TRUE.
2066                outflow_s = .FALSE.
2067             ENDIF
2068
2069             DEALLOCATE( f2_sub )
2070
2071          ELSE
2072
2073             CALL restrict_fast( f2, r )
2074
2075             ind_even_odd = even_odd_level(grid_level)  ! must be after restrict
2076
2077          ENDIF
2078
2079          p2 = 0.0_wp
2080
2081!
2082!--       Repeat the same procedure till the coarsest grid is reached
2083          CALL next_mg_level_fast( f2, p2, p3, r )
2084
2085       ENDIF
2086
2087!
2088!--    Now follows the prolongation
2089       IF ( grid_level >= 2 )  THEN
2090
2091!
2092!--       Prolongation of the new residual. The values are transferred
2093!--       from the coarse to the next finer grid.
2094          IF ( grid_level == mg_switch_to_pe0_level+1 )  THEN
2095
2096#if defined( __parallel )
2097!
2098!--          At this level, the new residual first has to be scattered from
2099!--          PE0 to the other PEs
2100             ALLOCATE( p2_sub(nzb:mg_loc_ind(5,myid)+1,             &
2101                       mg_loc_ind(3,myid)-1:mg_loc_ind(4,myid)+1,   &
2102                       mg_loc_ind(1,myid)-1:mg_loc_ind(2,myid)+1) )
2103
2104             CALL mg_scatter_fast( p2, p2_sub )
2105
2106!
2107!--          Therefore, indices of the previous level have to be changed to
2108!--          subdomain values in between (otherwise, the prolong routine would
2109!--          expect the gathered array)
2110             nxl_mg_save = nxl_mg(grid_level-1)
2111             nxr_mg_save = nxr_mg(grid_level-1)
2112             nys_mg_save = nys_mg(grid_level-1)
2113             nyn_mg_save = nyn_mg(grid_level-1)
2114             nzt_mg_save = nzt_mg(grid_level-1)
2115             nxl_mg(grid_level-1) = mg_loc_ind(1,myid)
2116             nxr_mg(grid_level-1) = mg_loc_ind(2,myid)
2117             nys_mg(grid_level-1) = mg_loc_ind(3,myid)
2118             nyn_mg(grid_level-1) = mg_loc_ind(4,myid)
2119             nzt_mg(grid_level-1) = mg_loc_ind(5,myid)
2120
2121!
2122!--          Set switch for routine exchange_horiz, that ghostpoint exchange
2123!--          has to be carried again out from now on
2124             mg_switch_to_pe0 = .FALSE.
2125
2126!
2127!--          For non-cyclic lateral boundary conditions, restore the
2128!--          in-/outflow conditions
2129             inflow_l  = .FALSE.;  inflow_r  = .FALSE.
2130             inflow_n  = .FALSE.;  inflow_s  = .FALSE.
2131             outflow_l = .FALSE.;  outflow_r = .FALSE.
2132             outflow_n = .FALSE.;  outflow_s = .FALSE.
2133
2134             IF ( pleft == MPI_PROC_NULL )  THEN
2135                IF ( bc_lr_dirrad )  THEN
2136                   inflow_l  = .TRUE.
2137                ELSEIF ( bc_lr_raddir )  THEN
2138                   outflow_l = .TRUE.
2139                ENDIF
2140             ENDIF
2141
2142             IF ( pright == MPI_PROC_NULL )  THEN
2143                IF ( bc_lr_dirrad )  THEN
2144                   outflow_r = .TRUE.
2145                ELSEIF ( bc_lr_raddir )  THEN
2146                   inflow_r  = .TRUE.
2147                ENDIF
2148             ENDIF
2149
2150             IF ( psouth == MPI_PROC_NULL )  THEN
2151                IF ( bc_ns_dirrad )  THEN
2152                   outflow_s = .TRUE.
2153                ELSEIF ( bc_ns_raddir )  THEN
2154                   inflow_s  = .TRUE.
2155                ENDIF
2156             ENDIF
2157
2158             IF ( pnorth == MPI_PROC_NULL )  THEN
2159                IF ( bc_ns_dirrad )  THEN
2160                   inflow_n  = .TRUE.
2161                ELSEIF ( bc_ns_raddir )  THEN
2162                   outflow_n = .TRUE.
2163                ENDIF
2164             ENDIF
2165
2166             CALL prolong_fast( p2_sub, p3 )
2167
2168!
2169!--          Restore the correct indices of the previous level
2170             nxl_mg(grid_level-1) = nxl_mg_save
2171             nxr_mg(grid_level-1) = nxr_mg_save
2172             nys_mg(grid_level-1) = nys_mg_save
2173             nyn_mg(grid_level-1) = nyn_mg_save
2174             nzt_mg(grid_level-1) = nzt_mg_save
2175
2176             DEALLOCATE( p2_sub )
2177#endif
2178
2179          ELSE
2180
2181             CALL prolong_fast( p2, p3 )
2182
2183          ENDIF
2184
2185!
2186!--       Computation of the new pressure correction. Therefore,
2187!--       values from prior grids are added up automatically stage by stage.
2188          DO  i = nxl_mg(grid_level)-1, nxr_mg(grid_level)+1
2189             DO  j = nys_mg(grid_level)-1, nyn_mg(grid_level)+1
2190                DO  k = nzb, nzt_mg(grid_level)+1
2191                   p_mg(k,j,i) = p_mg(k,j,i) + p3(k,j,i)
2192                ENDDO
2193             ENDDO
2194          ENDDO
2195
2196!
2197!--       Relaxation of the new solution
2198          CALL redblack_fast( f_mg, p_mg )
2199
2200       ENDIF
2201
2202
2203!
2204!--    The following few lines serve the steering of the multigrid scheme
2205       IF ( grid_level == maximum_grid_level )  THEN
2206
2207          GOTO 20
2208
2209       ELSEIF ( grid_level /= maximum_grid_level  .AND.  grid_level /= 1  .AND. &
2210                grid_level_count(grid_level) /= gamma_mg )  THEN
2211
2212          GOTO 10
2213
2214       ENDIF
2215
2216!
2217!--    Reset counter for the next call of poismg_fast
2218       grid_level_count(grid_level) = 0
2219
2220!
2221!--    Continue with the next finer level. nxl..nzt have to be
2222!--    set to the finer grid values, because these variables are needed for the
2223!--    exchange of ghost points in routine exchange_horiz
2224       grid_level = grid_level + 1
2225       ind_even_odd = even_odd_level(grid_level)
2226
2227       nxl = nxl_mg(grid_level)
2228       nxr = nxr_mg(grid_level)
2229       nys = nys_mg(grid_level)
2230       nyn = nyn_mg(grid_level)
2231       nzt = nzt_mg(grid_level)
2232
2233    20 CONTINUE
2234
2235    END SUBROUTINE next_mg_level_fast
2236
2237
2238
2239    SUBROUTINE init_even_odd_blocks
2240
2241!------------------------------------------------------------------------------!
2242! Description:
2243! ------------
2244! Initial settings for sorting k-dimension from sequential order (alternate
2245! even/odd) into blocks of even and odd or vice versa
2246!------------------------------------------------------------------------------!
2247
2248       USE arrays_3d,                                                          &
2249           ONLY:  f1_mg, f2_mg, f3_mg
2250
2251       USE control_parameters,                                                 &
2252           ONLY:  grid_level, masking_method, maximum_grid_level, topography
2253
2254       USE indices,                                                            &
2255           ONLY:  nzb, nzt, nzt_mg
2256
2257       USE indices,                                                            &
2258           ONLY:  flags, wall_flags_1, wall_flags_2, wall_flags_3,             &
2259                  wall_flags_4, wall_flags_5, wall_flags_6, wall_flags_7,      &
2260                  wall_flags_8, wall_flags_9, wall_flags_10, nxl_mg, nxr_mg,   &
2261                  nys_mg, nyn_mg, nzb, nzt_mg
2262
2263       IMPLICIT NONE
2264!
2265!--    Local variables
2266       INTEGER(iwp) ::  i     !:
2267       INTEGER(iwp) ::  l     !:
2268
2269       LOGICAL, SAVE ::  lfirst = .TRUE.
2270
2271
2272       IF ( .NOT. lfirst )  RETURN
2273
2274       ALLOCATE( even_odd_level(maximum_grid_level) )
2275
2276       ALLOCATE( f1_mg_b(nzb:nzt+1,maximum_grid_level),                        &
2277                 f2_mg_b(nzb:nzt+1,maximum_grid_level),                        &
2278                 f3_mg_b(nzb:nzt+1,maximum_grid_level) )
2279
2280!
2281!--    Set border index between the even and odd block
2282       DO  i = maximum_grid_level, 1, -1
2283          even_odd_level(i) = nzt_mg(i) / 2
2284       ENDDO
2285
2286!
2287!--    Sort grid coefficients used in red/black scheme and for calculating the
2288!--    residual to block (even/odd) structure
2289       DO  l = maximum_grid_level, 1 , -1
2290          CALL sort_k_to_even_odd_blocks( f1_mg(nzb+1:nzt_mg(grid_level),l),   &
2291                                          f1_mg_b(nzb:nzt_mg(grid_level)+1,l), &
2292                                          l )
2293          CALL sort_k_to_even_odd_blocks( f2_mg(nzb+1:nzt_mg(grid_level),l),   &
2294                                          f2_mg_b(nzb:nzt_mg(grid_level)+1,l), &
2295                                          l )
2296          CALL sort_k_to_even_odd_blocks( f3_mg(nzb+1:nzt_mg(grid_level),l),   &
2297                                          f3_mg_b(nzb:nzt_mg(grid_level)+1,l), &
2298                                          l )
2299       ENDDO
2300
2301!
2302!--    Sort flags array for all levels to block (even/odd) structure
2303       IF ( topography /= 'flat'  .AND.  .NOT. masking_method )  THEN
2304
2305          DO  l = maximum_grid_level, 1 , -1
2306!
2307!--          Assign the flag level to be calculated
2308             SELECT CASE ( l )
2309                CASE ( 1 )
2310                   flags => wall_flags_1
2311                CASE ( 2 )
2312                   flags => wall_flags_2
2313                CASE ( 3 )
2314                   flags => wall_flags_3
2315                CASE ( 4 )
2316                   flags => wall_flags_4
2317                CASE ( 5 )
2318                   flags => wall_flags_5
2319                CASE ( 6 )
2320                   flags => wall_flags_6
2321                CASE ( 7 )
2322                   flags => wall_flags_7
2323                CASE ( 8 )
2324                   flags => wall_flags_8
2325                CASE ( 9 )
2326                   flags => wall_flags_9
2327                CASE ( 10 )
2328                   flags => wall_flags_10
2329             END SELECT
2330
2331             CALL sort_k_to_even_odd_blocks( flags, l )
2332
2333          ENDDO
2334
2335       ENDIF
2336
2337       lfirst = .FALSE.
2338
2339     END SUBROUTINE init_even_odd_blocks
2340
2341
2342
2343     SUBROUTINE special_exchange_horiz ( p_mg, color )
2344
2345!------------------------------------------------------------------------------!
2346! Description:
2347! ------------
2348! Special exchange_horiz subroutine for use in redblack. Transfers only
2349! "red" or "black" data points.
2350!------------------------------------------------------------------------------!
2351
2352       USE control_parameters,                                                 &
2353           ONLY:  bc_lr_cyc, bc_ns_cyc, grid_level, ibc_p_b, ibc_p_t,          &
2354                  inflow_l, inflow_n, inflow_r, inflow_s, maximum_grid_level,  &
2355                  outflow_l, outflow_n, outflow_r, outflow_s,                  &
2356                  synchronous_exchange
2357
2358       USE indices,                                                            &
2359           ONLY:  mg_loc_ind, nxl, nxl_mg, nxr, nxr_mg, nys, nys_mg, nyn,      &
2360                  nyn_mg, nzb, nzt, nzt_mg
2361
2362       IMPLICIT NONE
2363
2364       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level)+1,                           &
2365                           nys_mg(grid_level)-1:nyn_mg(grid_level)+1,          &
2366                           nxl_mg(grid_level)-1:nxr_mg(grid_level)+1) ::  p_mg  !:
2367
2368       INTEGER(iwp), intent(IN) ::  color  !:
2369!
2370!--    Local variables
2371       INTEGER(iwp) ::  i,i1,i2  !:
2372       INTEGER(iwp) ::  j,j1,j2  !:
2373       INTEGER(iwp) ::  k        !:
2374       INTEGER(iwp) ::  l        !:
2375       INTEGER(iwp) ::  jys      !:
2376       INTEGER(iwp) ::  jyn      !:
2377       INTEGER(iwp) ::  ixl      !:
2378       INTEGER(iwp) ::  ixr      !:
2379       logical      ::  sendrecv_in_background_save  !:
2380       logical      ::  synchronous_exchange_save    !:
2381
2382       REAL(wp), DIMENSION(nzb:nzt_mg(grid_level-1)+1,                         &
2383                           nys_mg(grid_level-1)-1:nyn_mg(grid_level-1)+1,      &
2384                           nxl_mg(grid_level-1)-1:nxr_mg(grid_level-1)+1) ::  temp  !:
2385
2386#if defined ( __parallel )
2387       sendrecv_in_background_save = sendrecv_in_background
2388       sendrecv_in_background      = .FALSE.
2389       synchronous_exchange_save   = synchronous_exchange
2390       synchronous_exchange        = .FALSE.
2391
2392       l = grid_level
2393
2394       ind_even_odd = even_odd_level(grid_level)
2395
2396       jys = nys_mg(grid_level-1)
2397       jyn = nyn_mg(grid_level-1)
2398       ixl = nxl_mg(grid_level-1)
2399       ixr = nxr_mg(grid_level-1)
2400
2401!
2402!--    Restricted transfer only on finer levels with enough data
2403       IF ( ngp_xz(grid_level) >= 900 .OR. ngp_yz(grid_level) >= 900 )  THEN
2404!
2405!--       Handling the even k Values
2406!--       Collecting data for the north - south exchange
2407!--       Since only every second value has to be transfered, data are stored
2408!--       on the next coarser grid level, because the arrays on that level
2409!--       have just the required size
2410          i1 = nxl_mg(grid_level-1)
2411          i2 = nxl_mg(grid_level-1)
2412
2413          DO  i = nxl_mg(l), nxr_mg(l), 2
2414             DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2415
2416                IF ( j == nys_mg(l) )  THEN
2417                   !DIR$ IVDEP
2418                   DO  k = ind_even_odd+1, nzt_mg(l)
2419                      temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
2420                   ENDDO
2421                   i1 = i1 + 1
2422
2423                ENDIF
2424
2425                IF ( j == nyn_mg(l) )  THEN
2426                   !DIR$ IVDEP
2427                   DO  k = ind_even_odd+1, nzt_mg(l)
2428                      temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
2429                   ENDDO
2430                   i2 = i2 + 1
2431
2432                ENDIF
2433
2434             ENDDO
2435          ENDDO
2436
2437          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2438             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2439
2440                IF ( j == nys_mg(l) )  THEN
2441                   !DIR$ IVDEP
2442                   DO  k = ind_even_odd+1, nzt_mg(l)
2443                      temp(k-ind_even_odd,jys,i1) = p_mg(k,j,i)
2444                   ENDDO
2445                   i1 = i1 + 1
2446
2447                ENDIF
2448
2449                IF ( j == nyn_mg(l) )  THEN
2450                   !DIR$ IVDEP
2451                   DO  k = ind_even_odd+1, nzt_mg(l)
2452                      temp(k-ind_even_odd,jyn,i2) = p_mg(k,j,i)
2453                   ENDDO
2454                   i2 = i2 + 1
2455
2456                ENDIF
2457
2458             ENDDO
2459          ENDDO
2460
2461          grid_level = grid_level-1
2462
2463          nxl = nxl_mg(grid_level)
2464          nys = nys_mg(grid_level)
2465          nxr = nxr_mg(grid_level)
2466          nyn = nyn_mg(grid_level)
2467          nzt = nzt_mg(grid_level)
2468
2469          send_receive = 'ns'
2470          CALL exchange_horiz( temp, 1 )
2471
2472          grid_level = grid_level+1
2473
2474          i1 = nxl_mg(grid_level-1)
2475          i2 = nxl_mg(grid_level-1)
2476
2477          DO  i = nxl_mg(l), nxr_mg(l), 2
2478             DO  j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2479
2480                IF ( j == nys_mg(l) )  THEN
2481                   !DIR$ IVDEP
2482                   DO  k = ind_even_odd+1, nzt_mg(l)
2483                      p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
2484                   ENDDO
2485                   i1 = i1 + 1
2486
2487                ENDIF
2488
2489                IF ( j == nyn_mg(l) )  THEN
2490                   !DIR$ IVDEP
2491                   DO  k = ind_even_odd+1, nzt_mg(l)
2492                      p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
2493                   ENDDO
2494                   i2 = i2 + 1
2495
2496                ENDIF
2497
2498             ENDDO
2499          ENDDO
2500
2501          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2502             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2503
2504                IF ( j == nys_mg(l) )  THEN
2505                   !DIR$ IVDEP
2506                   DO  k = ind_even_odd+1, nzt_mg(l)
2507                      p_mg(k,nyn_mg(l)+1,i) = temp(k-ind_even_odd,jyn+1,i1)
2508                   ENDDO
2509                   i1 = i1 + 1
2510
2511                ENDIF
2512
2513                IF ( j == nyn_mg(l) )  THEN
2514                   !DIR$ IVDEP
2515                   DO  k = ind_even_odd+1, nzt_mg(l)
2516                      p_mg(k,nys_mg(l)-1,i) = temp(k-ind_even_odd,jys-1,i2)
2517                   ENDDO
2518                   i2 = i2 + 1
2519
2520                ENDIF
2521
2522             ENDDO
2523          ENDDO
2524
2525!
2526!--       Collecting data for the left - right exchange
2527!--       Since only every second value has to be transfered, data are stored
2528!--       on the next coarser grid level, because the arrays on that level
2529!--       have just the required size
2530          j1 = nys_mg(grid_level-1)
2531          j2 = nys_mg(grid_level-1)
2532
2533          DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2534             DO  i = nxl_mg(l), nxr_mg(l), 2
2535
2536                IF ( i == nxl_mg(l) )  THEN
2537                   !DIR$ IVDEP
2538                   DO  k = ind_even_odd+1, nzt_mg(l)
2539                      temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
2540                   ENDDO
2541                   j1 = j1 + 1
2542
2543                ENDIF
2544
2545                IF ( i == nxr_mg(l) )  THEN
2546                   !DIR$ IVDEP
2547                   DO  k = ind_even_odd+1, nzt_mg(l)
2548                      temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
2549                   ENDDO
2550                   j2 = j2 + 1
2551
2552                ENDIF
2553
2554             ENDDO
2555          ENDDO
2556
2557          DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2558             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2559
2560                IF ( i == nxl_mg(l) )  THEN
2561                   !DIR$ IVDEP
2562                   DO  k = ind_even_odd+1, nzt_mg(l)
2563                      temp(k-ind_even_odd,j1,ixl) = p_mg(k,j,i)
2564                   ENDDO
2565                   j1 = j1 + 1
2566
2567                ENDIF
2568
2569                IF ( i == nxr_mg(l) )  THEN
2570                   !DIR$ IVDEP
2571                   DO  k = ind_even_odd+1, nzt_mg(l)
2572                      temp(k-ind_even_odd,j2,ixr) = p_mg(k,j,i)
2573                   ENDDO
2574                   j2 = j2 + 1
2575
2576                ENDIF
2577
2578             ENDDO
2579          ENDDO
2580
2581          grid_level = grid_level-1
2582          send_receive = 'lr'
2583
2584          CALL exchange_horiz( temp, 1 )
2585
2586          grid_level = grid_level+1
2587
2588          j1 = nys_mg(grid_level-1)
2589          j2 = nys_mg(grid_level-1)
2590
2591          DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2592             DO  i = nxl_mg(l), nxr_mg(l), 2
2593
2594                IF ( i == nxl_mg(l) )  THEN
2595                   !DIR$ IVDEP
2596                   DO  k = ind_even_odd+1, nzt_mg(l)
2597                      p_mg(k,j,nxr_mg(l)+1)  = temp(k-ind_even_odd,j1,ixr+1)
2598                   ENDDO
2599                   j1 = j1 + 1
2600
2601                ENDIF
2602
2603                IF ( i == nxr_mg(l) )  THEN
2604                   !DIR$ IVDEP
2605                   DO  k = ind_even_odd+1, nzt_mg(l)
2606                      p_mg(k,j,nxl_mg(l)-1)  = temp(k-ind_even_odd,j2,ixl-1)
2607                   ENDDO
2608                   j2 = j2 + 1
2609
2610                ENDIF
2611
2612             ENDDO
2613          ENDDO
2614
2615          DO j = nys_mg(l) + (color-1), nyn_mg(l), 2
2616             DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2617
2618                IF ( i == nxl_mg(l) )  THEN
2619                   !DIR$ IVDEP
2620                   DO  k = ind_even_odd+1, nzt_mg(l)
2621                      p_mg(k,j,nxr_mg(l)+1)  = temp(k-ind_even_odd,j1,ixr+1)
2622                   ENDDO
2623                   j1 = j1 + 1
2624
2625                ENDIF
2626
2627                IF ( i == nxr_mg(l) )  THEN
2628                   !DIR$ IVDEP
2629                   DO  k = ind_even_odd+1, nzt_mg(l)
2630                      p_mg(k,j,nxl_mg(l)-1)  = temp(k-ind_even_odd,j2,ixl-1)
2631                   ENDDO
2632                   j2 = j2 + 1
2633
2634                ENDIF
2635
2636             ENDDO
2637          ENDDO
2638
2639!
2640!--       Now handling the even k values
2641!--       Collecting data for the north - south exchange
2642!--       Since only every second value has to be transfered, data are stored
2643!--       on the next coarser grid level, because the arrays on that level
2644!--       have just the required size
2645          i1 = nxl_mg(grid_level-1)
2646          i2 = nxl_mg(grid_level-1)
2647
2648          DO  i = nxl_mg(l), nxr_mg(l), 2
2649             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2650
2651                IF ( j == nys_mg(l) )  THEN
2652                   !DIR$ IVDEP
2653                   DO  k = nzb+1, ind_even_odd
2654                      temp(k,jys,i1) = p_mg(k,j,i)
2655                   ENDDO
2656                   i1 = i1 + 1
2657
2658                ENDIF
2659
2660                IF ( j == nyn_mg(l) )  THEN
2661                   !DIR$ IVDEP
2662                   DO  k = nzb+1, ind_even_odd
2663                      temp(k,jyn,i2) = p_mg(k,j,i)
2664                   ENDDO
2665                   i2 = i2 + 1
2666
2667                ENDIF
2668
2669             ENDDO
2670          ENDDO
2671
2672          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2673             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2674
2675                IF ( j == nys_mg(l) )  THEN
2676                   !DIR$ IVDEP
2677                   DO  k = nzb+1, ind_even_odd
2678                      temp(k,jys,i1) = p_mg(k,j,i)
2679                   ENDDO
2680                   i1 = i1 + 1
2681
2682                ENDIF
2683
2684                IF ( j == nyn_mg(l) )  THEN
2685                   !DIR$ IVDEP
2686                   DO  k = nzb+1, ind_even_odd
2687                      temp(k,jyn,i2) = p_mg(k,j,i)
2688                   ENDDO
2689                   i2 = i2 + 1
2690
2691                ENDIF
2692
2693             ENDDO
2694          ENDDO
2695
2696          grid_level = grid_level-1
2697
2698          send_receive = 'ns'
2699          CALL exchange_horiz( temp, 1 )
2700
2701          grid_level = grid_level+1
2702
2703          i1 = nxl_mg(grid_level-1)
2704          i2 = nxl_mg(grid_level-1)
2705
2706          DO  i = nxl_mg(l), nxr_mg(l), 2
2707             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2708
2709                IF ( j == nys_mg(l) )  THEN
2710                   !DIR$ IVDEP
2711                   DO  k = nzb+1, ind_even_odd
2712                      p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
2713                   ENDDO
2714                   i1 = i1 + 1
2715
2716                ENDIF
2717
2718                IF ( j == nyn_mg(l) )  THEN
2719                   !DIR$ IVDEP
2720                   DO  k = nzb+1, ind_even_odd
2721                      p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
2722                   ENDDO
2723                   i2 = i2 + 1
2724
2725                ENDIF
2726
2727             ENDDO
2728          ENDDO
2729
2730          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2731             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2732
2733                IF ( j == nys_mg(l) )  THEN
2734                   !DIR$ IVDEP
2735                   DO  k = nzb+1, ind_even_odd
2736                      p_mg(k,nyn_mg(l)+1,i) = temp(k,jyn+1,i1)
2737                   ENDDO
2738                   i1 = i1 + 1
2739
2740                ENDIF
2741
2742                IF ( j == nyn_mg(l) )  THEN
2743                   !DIR$ IVDEP
2744                   DO  k = nzb+1, ind_even_odd
2745                      p_mg(k,nys_mg(l)-1,i) = temp(k,jys-1,i2)
2746                   ENDDO
2747                   i2 = i2 + 1
2748
2749                ENDIF
2750
2751             ENDDO
2752          ENDDO
2753
2754          j1 = nys_mg(grid_level-1)
2755          j2 = nys_mg(grid_level-1)
2756
2757          DO  i = nxl_mg(l), nxr_mg(l), 2
2758             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2759
2760                IF ( i == nxl_mg(l) )  THEN
2761                   !DIR$ IVDEP
2762                   DO  k = nzb+1, ind_even_odd
2763                      temp(k,j1,ixl) = p_mg(k,j,i)
2764                   ENDDO
2765                   j1 = j1 + 1
2766
2767                ENDIF
2768
2769                IF ( i == nxr_mg(l) )  THEN
2770                   !DIR$ IVDEP
2771                   DO  k = nzb+1, ind_even_odd
2772                      temp(k,j2,ixr) = p_mg(k,j,i)
2773                   ENDDO
2774                   j2 = j2 + 1
2775
2776                ENDIF
2777
2778             ENDDO
2779          ENDDO
2780
2781          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2782             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2783
2784                IF ( i == nxl_mg(l) )  THEN
2785                   !DIR$ IVDEP
2786                   DO  k = nzb+1, ind_even_odd
2787                      temp(k,j1,ixl) = p_mg(k,j,i)
2788                   ENDDO
2789                   j1 = j1 + 1
2790
2791                ENDIF
2792
2793                IF ( i == nxr_mg(l) )  THEN
2794                   !DIR$ IVDEP
2795                   DO  k = nzb+1, ind_even_odd
2796                      temp(k,j2,ixr) = p_mg(k,j,i)
2797                   ENDDO
2798                   j2 = j2 + 1
2799
2800                ENDIF
2801
2802             ENDDO
2803          ENDDO
2804
2805          grid_level = grid_level-1
2806
2807          send_receive = 'lr'
2808          CALL exchange_horiz( temp, 1 )
2809
2810          grid_level = grid_level+1
2811
2812          nxl = nxl_mg(grid_level)
2813          nys = nys_mg(grid_level)
2814          nxr = nxr_mg(grid_level)
2815          nyn = nyn_mg(grid_level)
2816          nzt = nzt_mg(grid_level)
2817
2818          j1 = nys_mg(grid_level-1)
2819          j2 = nys_mg(grid_level-1)
2820
2821          DO  i = nxl_mg(l), nxr_mg(l), 2
2822             DO  j = nys_mg(l) + (color-1), nyn_mg(l), 2
2823
2824                IF ( i == nxl_mg(l) )  THEN
2825                   !DIR$ IVDEP
2826                   DO  k = nzb+1, ind_even_odd
2827                      p_mg(k,j,nxr_mg(l)+1)  = temp(k,j1,ixr+1)
2828                   ENDDO
2829                   j1 = j1 + 1
2830
2831                ENDIF
2832
2833                IF ( i == nxr_mg(l) )  THEN
2834                   !DIR$ IVDEP
2835                   DO  k = nzb+1, ind_even_odd
2836                      p_mg(k,j,nxl_mg(l)-1)  = temp(k,j2,ixl-1)
2837                   ENDDO
2838                   j2 = j2 + 1
2839
2840                ENDIF
2841
2842             ENDDO
2843          ENDDO
2844
2845          DO  i = nxl_mg(l)+1, nxr_mg(l), 2
2846             DO j = nys_mg(l) + 2 - color, nyn_mg(l), 2
2847
2848                IF ( i == nxl_mg(l) )  THEN
2849                   !DIR$ IVDEP
2850                   DO  k = nzb+1, ind_even_odd
2851                      p_mg(k,j,nxr_mg(l)+1)  = temp(k,j1,ixr+1)
2852                   ENDDO
2853                   j1 = j1 + 1
2854
2855                ENDIF
2856
2857                IF ( i == nxr_mg(l) )  THEN
2858                   !DIR$ IVDEP
2859                   DO  k = nzb+1, ind_even_odd
2860                      p_mg(k,j,nxl_mg(l)-1)  = temp(k,j2,ixl-1)
2861                   ENDDO
2862                   j2 = j2 + 1
2863
2864                ENDIF
2865
2866             ENDDO
2867          ENDDO
2868
2869       ELSE
2870
2871!
2872!--       Standard horizontal ghost boundary exchange for small coarse grid
2873!--       levels, where the transfer time is latency bound
2874          CALL exchange_horiz( p_mg, 1 )
2875
2876       ENDIF
2877
2878!
2879!--    Reset values to default PALM setup
2880       sendrecv_in_background = sendrecv_in_background_save
2881       synchronous_exchange   = synchronous_exchange_save
2882       send_receive = 'al'
2883#else
2884
2885!
2886!--    Standard horizontal ghost boundary exchange for small coarse grid
2887!--    levels, where the transfer time is latency bound
2888       CALL exchange_horiz( p_mg, 1 )
2889#endif
2890
2891
2892
2893    END SUBROUTINE special_exchange_horiz
2894
2895 END MODULE poismg_mod
Note: See TracBrowser for help on using the repository browser.