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

Last change on this file since 1807 was 1763, checked in by hellstea, 8 years ago

last commit documented

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