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

Last change on this file since 1850 was 1850, checked in by maronga, 6 years ago

added _mod string to several filenames to meet the naming convection for modules

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