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

Last change on this file since 1762 was 1762, checked in by hellstea, 6 years ago

Introduction of nested domain system

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