source: palm/trunk/SOURCE/poismg_fast_mod.f90 @ 1856

Last change on this file since 1856 was 1851, checked in by maronga, 8 years ago

last commit documented

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