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

Last change on this file since 1683 was 1683, checked in by knoop, 6 years ago

last commit documented

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